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1  Executive  Summary 


This  final  report  titled  “STORESIM:  An  Integrated  System  for  Multi-Body  CFD  Simulation 
Using  Unstructured,  Adaptive  Grids”  documents  the  effort  conducted  at  the  Computational 
Mechanics  Company,  Inc.,  (COMCO)  under  contract  No.  F33615-90-C-3001  with  the  Air 
Force  Flight  Dynamics  Directorate,  Wright  Laboratory,  OH  45433-7562.  The  work  spans  from 
September  1990  to  August  1996.  The  Air  Force  Monitors  for  the  project  were  F.  Witzeman 
and  J.  Nelson.  The  work  was  conducted  by  COMCO  in  Austin,  Texas.  Key  investigators  were 
Dr.  Stephen  R.  Kennon,  Dr.  Chittur  S.  Venkatasubban,  Mr.  Jim  M.  Meyering,  Mr.  Manas 
K.  Deb,  and  Dr.  Mahender  P.  Reddy. 

The  aerodynamic  simulation  of  store  separation,  missile  launching,  and  crew  escape  is  a 
mission  critical  problem  class  which  only  recently  has  become  amenable  to  computational 
fluid  dynamics  (CFD)  modeling  techniques  coupled  with  high  performance  computing.  These 
problems  fall  under  the  general  classification  of  Relative  Body  Motion  Aerodynamics  (RBMA), 
characterized  by  multiple  moving  bodies  immersed  in  a  viscous  fluid.  CFD  provides  a  unique, 
essential  tool  for  analyzing  these  difficult  problems  allowing  the  study  of  arbitrary  configura¬ 
tions  without  the  expense  of  model  building,  wind-tunnel  testing,  and/or  flight  testing. 

Numerical  analysis  of  RBMA  problems  involves  grid  generation,  grid  movement/  restruc¬ 
turing,  adaptive  methods,  flow  simulation,  high  temperature  gasdynamics,  turbulence,  rigid 
body  dynamics,  and  interactive  pre  and  post-processing.  The  primary  goal  of  this  effort  was  to 
couple  these  disciplines  into  a  unified  computational  environment  for  analyzing  dynamic  sep¬ 
aration  of  weapons,  stages,  and  crew  escape  vehicles  from  supersonic  and  hypersonic  parent 
vehicles.  The  major  outcome  of  this  research  and  development  effort  is  a  software  pack¬ 
age  called  STORESIM.  The  STORESIM  package  is  a  tightly  integrated  group  of  modules 
(TETMESH,  STOREFLOW,  STOREVIS)  that  performs  grid  generation,  static  flow  simu¬ 
lation,  six  degree-of-freedom  rigid  body  dynamics,  dynamic  flow  analysis,  body  movement, 
grid-node  movement,  grid  restructuring,  and  time-dependent  results  visualization. 

Key  accomplishments  of  this  effort  are: 

1.  Mesh  Generation:  A  grid  generator  is  essentially  useless  without  a  powerful  technique  for 
defining  the  geometry  to  be  discretized.  To  this  end  we  have  developed  an  unstructured 
mesh  generation  module,  TETMESH,  that  is  capable  of  reading  the  geometry  in  a  very 
simple  network/trimmed  surface  format.  Salient  features  of  this  module  include: 

•  Triangulation:  Surface  triangulation  is  based  on  the  Dulanay  criterion.  The  surface 
triangulator  is  capable  of  handling  complex  curved  surfaces  with  cusped  edges. 

•  Volume  Gridding:  The  volume  mesher  generates  tetrahedral  elements  with  an  op¬ 
tion  of  generating  prismatic  elements  extruded  from  the  surfaces.  The  tetrahedral- 
ization  is  also  based  on  Dulanay  criterion  with  proper  face  enforcement  to  fill  the 
voids. 

•  Zero-thickness  Surface:  Thin  fins  protruding  from  a  missile  can  be  approximated 
as  single  surfaces  during  data  preparation.  TETMESH  automatically  generates  a 
matching  surface  to  create  a  barrier  and  prevent  fluid  movement  across  this  surface. 
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•  Adaptivity:  The  meshing  module  is  capable  of  refining/unrefining  the  mesh.  Re¬ 
gions  of  maximum  solution  gradients  are  automatically  refined  to  better  capture 
solution  in  critical  regions. 

2.  Flow  Solver:  The  STORESIM  package  has  a  compressible  fluid  flow  solver  module  called 
STOREFLOW.  Key  features  of  the  STOREFLOW  module  are: 

•  An  option  for  ALE  (Arbitrary  Lagrangian  Eulerian )  computations.  This  essentially 
makes  it  possible  to  use  moving  meshes  to  model  multiple  moving  bodies.  This 
new  ALE  formulation  is  extremely  robust,  and  can  even  produce  accurate  results 
for  bodies  launched  at  high  speed  into  ambient  air. 

•  A  very  sophisticated  far  stream  boundary  condition  designed  to  prevent  spurious 
energy  build  up  and  consequent  reflections  at  transonic  speeds.  This  approach  also 
includes  modifications  to  account  for  moving  mesh  boundaries. 

•  Euler  and  Navier-Stokes  versions.  Turbulence  modeling  based  on  a  simple  zero 
equation  model  has  been  tested  for  two-dimensional  flows.  The  extension  to  three- 
dimensional  turbulence  modeling,  although  not  implements  should  be  straightfor¬ 
ward. 

•  Real  gas  capability,  based  on  a  5  species  (Nitrogen,  Oxygen,  Atomic  Nitrogen, 
Atomic  Oxygen  and  Nitric  Oxide)  model  for  dissociated  air. 

•  Multiple  element  types  (tetrahedra,  Hexahedra  or  Prisms).  Additional  element 
types  can  be  easily  incorporated  into  the  code. 

3.  Six-DOF  Solver:  STOREFLOW  is  fully  integrated  with  a  six  of  freedom  (6  dof)  solver 
for  rigid  body  dynamics,  and  a  meshmover,  that  moves  the  mesh  in  response  to  body 
motions. 

4.  Node  Movement:  An  innovative  grid  movement /restructuring  method  based  on  the  grid- 
node-flow  analogy  and  generalized  edge  swapping  in  three-dimensions.  The  node  move¬ 
ment/restructuring  is  required  when  modeling  body  movements  in  the  flow  field. 

5.  Visualization:  The  visualization  package,  STOREVIS,  is  capable  of  post-processing 
structured/unstructured  grids  and  producing  color  plots  of  slices,  isosurfaces,  velocity 
vectors,  and  particle  traces. 

STORESIM  has  been  used  to  solve  several  complex  problems.  In  every  case,  all  relevant 
modules  of  the  STORESIM  package  are  used.  During  the  development  and  final  validation 
process,  each  module  has  been  tested  separately  and  as  a  complete  package  for  accuracy  and 
robustness.  The  final  code  has  been  transferred  to  Air  Force  Flight  Dynamics  Directorate. 
This  report  discusses  these  efforts  and  describes  the  STORESIM  module  in  detail. 
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2  Introduction 


The  aerodynamic  simulation  of  store  separation,  missile  launching,  and  crew  escape  is  a  mis¬ 
sion  critical  problem  class  which  only  recently  has  become  amenable  to  computational  fluid 
dynamics  (CFD)  modeling  techniques  coupled  with  high  performance  computing.  These  prob¬ 
lems  fall  under  the  general  classification  of  Relative  Body  Motion  Aerodynamics  (RBMA), 
characterized  by  multiple  moving  bodies  immersed  in  a  viscous  fluid.  CFD  provides  a  unique, 
essential  tool  for  analyzing  these  difficult  problems  allowing  the  study  of  arbitrary  configura¬ 
tions  without  the  expense  of  model  building,  wind-tunnel  testing,  and/or  flight  testing. 

When  multiple  bodies  that  are  moving  relative  to  each  other  are  immersed  in  a  high-speed 
fluid  stream,  many  complex  flow  phenomena  exist  that  must  be  carefully  modeled.  In  many 
situations,  the  bodies  start  out  or  become  close  together  producing  strong  shocks  in  choking 
flow  associated  with  shock/viscous  interactions.  These  choking  shocks  can  produce  adverse 
aerothermodynamic  affects  that  when  coupled  with  the  overall  proximity  to  other  bodies,  can 
result  in  undesirable  flight  trajectories  of  missiles,  stores,  and/or  crew  escape  vehicles.  In 
some  situations,  these  adverse  aerodynamic  effects  can  be  severe  enough  to  produce  contact 
between  the  separating  body  and  the  parent  aircraft.  The  study  of  these  effects  is  even  more 
important  for  high-speed  supersonic  and  hypersonic  aircraft  of  the  future.  These  aircraft  must 
be  able  to  separate  safely  from  their  stores  and/or  missiles  at  high  speeds. 

In  addition  to  the  proximity  problems  associated  with  one  or  more  bodies,  the  effects 
of  unsteady  flows  are  important  in  RBMA.  Shocks,  shear  layers  and  vortices  appear  and 
disappear  within  a  RBMA  simulation.  The  seldom  seen  strong  oblique  shock  solution  of  the 
gasdynamic  equations  can  appear  due  to  the  motion  of  a  body  relative  to  the  free-stream  flow. 
To  accurately  model  these  diverse  affects  requires  the  use  of  modern  adaptive  grid  technologies 
that  are  able  to  adapt  the  grid  to  regions  of  complex  unsteady  flow  features. 

Another  phenomena  which  will  encountered  by  future  aircraft  is  hypersonic  flow-fields 
with  associated  effects.  These  flows  are  characterized  by  high  temperature  effects  that  must 
be  accurately  modeled  and  present  still  another  challenge  to  modern  CFD.  Finally,  the  effects 
of  turbulence  must  be  accounted  for  by  appropriate  models. 

2.1  Objectives 

The  US  Air  Force  has  a  great  need  for  modern  CFD  simulation  software  that  will  quickly 
and  accurately  predict  the  behavior  of  new  and/or  modified  weapons  systems.  This  need  has 
arisen  due  to  the  high  cost  of  wind-tunnel-  and  flight-testing  of  new  weapons  configurations  on 
current  and  future  aircraft.  Numerical  analysis  of  RBMA  problems  involves  grid  generation, 
grid  movement/  restructuring,  adaptive  methods,  flow  solution,  high  temperature  gasd3mam- 
ics,  turbulence,  rigid  body  dynamics,  and  interactive  pre  and  post-processing.  A  successful 
simulation  environment  requires  an  integrated  package  that  is  easy  to  use,  can  accept  standard 
CAD  geometries  from  USAF  databases,  and  is  automatic  and  requires  minimal  user  input  or 
steering. 

The  main  objective  of  this  effort  was  to  develop  one  such  state-of-the-art  software  package 
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to  accurately  and  efficiently  model  these  diverse  phenomena.  To  this  end  we  have  developed 
STORESIM.  As  mentioned  earlier,  STORESIM  is  an  integrated  package  of  modules  that 
performs  grid  generation,  static  flow  solution,  six  degree-of-freedora  rigid  body  dynamics, 
dynamic  flow  solution,  body  movement,  grid  node  movement,  grid  restructuring,  and  time- 
dependent  results  visualization.  Although  the  package  is  tightly  integrated,  and  presents 
a  uniform  view  of  the  simulation  to  the  user,  the  software  architecture  is  based  on  a  set 
of  modules  that  can  be  used  as  stand-alone  programs.  This  report  describes  each  of  these 
modules  in  detail.  In  the  next  Section  we  describe  the  mesh  generation  procedure  using 
the  meshing  module  called  TETMESH.  In  Sections  3  we  describe  the  flow  solver  and  the 
modifications  related  to  implementation  of  the  ALE  scheme.  In  Section  4  we  explain  the 
six-dof  rigid-body-dynamics  solver  and  the  mesh  movement  algorithm.  The  flow  visualization 
package  (STOREVIS)  is  explained  in  Section  5.  Finally  in  Sections  6  and  7  we  present  the 
validation  cases  and  demonstration  cases  showing  vaious  usages  of  STORESIM  in  typical 
problems  encountered  in  real-life  applications. 


4 


3  Mesh  Generation 


The  efficient  generation  of  superior  quality  computational  grids  for  arbitrary  configurations  has 
long  been  a  pacing  item  in  CFD.  The  advent  of  unstructured  tetrahedral  grids  has  made  the 
grid  generation  task  an  order  of  magnitude  easier  than  either  structured  or  block-structured 
approaches  due  to  the  inherent  generality  of  unstructured  grids.  In  the  following  we  explain 
the  capabilities  of  TETMESH,  the  STORESIM  mesh  generation  module.  Before  discussing 
TETMESH,  we  digress  to  discuss  the  general  three-dimensional  tetrahedralization  problem. 

3.1  The  3D  Tetrahedralization  Problem 

The  problem  can  be  stated  as  follows: 

“Given  a  set  of  triangles  whose  union  forms  the  boundary  for  a  closed  region  in 
3D,  fill  the  region  with  tetrahedra  such  that  the  surface  triangles  appear  as  faces  of 
the  tetrahedralization.  ” 

At  first  glance  (and  first  implementation  as  well),  one  might  think  that  this  problem  is 
trivial.  However,  the  basic  fact  is  that  unstructured  mesh  generation  in  three  dimensions  is 
orders  of  magnitude  more  difficult  than  the  analogous  problem  in  two  dimensions.  This  is  due 
to  basic  theorems  that  state  that  a  given  polyhedron  formed  from  triangular  facets  does  not 
always  have  a  valid  tetrahedralization  [19]. 

The  basic  methods  for  solving  the  above  problem  are  Advancing  Front  and  Delaunay  (Oct- 
tree  solves  a  different  problem,  since  the  method  creates  the  surface  triangles/nodes  as  part  of 
the  process).  The  Advancing  Front  method  is  appealing  since  one  can  practically  guarantee 
that  one  can  march  the  mesh  away  from  the  surface  for  at  least  a  few  layers  and,  therefore,  a 
priori  enforce  that  the  surface  triangles  appear.  However,  one  has  simply  pushed  the  problem 
away  from  the  boundary  into  the  interior  by  forming  a  (possibly  more  complicated)  new  cavity 
that  must  be  filled  with  tetrahedra.  Eventually,  as  borne  out  by  practical  experience  with  these 
methods,  a  non-tetrahedralizable  cavity  is  often  produced,  and  the  standard  method  will  fail 
since  it  will  not  be  able  to  fill  the  last  hole  in  the  mesh  (note:  various  schemes  are  used  to 
try  to  overcome  this,  including  backing  up  a  few  layers,  etc.).  The  Delaunay  method  is  no 
different  in  that  it  suffers  from  the  problem  right  at  the  domain  boundaries,  except  for  the 
few  simple  cases  where  the  surface  triangulation  is  Delaunay  compatible  (see  below)  with  the 
volume  mesh.  The  conclusion  is  that  for  either  approach,  one  will  need  (eventually)  a  Face 
Enforcer  that  can  enforce  the  required  cavity  faces  100%  of  the  time.  We  have  implemented 
one  such  face  enforcer  and  it  is  described  below. 

The  module  in  STORESIM  for  grid  generation,  TETMESH,  implements  a  constrained 
Delaunay  procedure  for  surface  triangulation  and  interior  tetrahedralization.  TETMESH  also 
generates  semi-structured  prismatic  grids  by  extruding  surface  triangles  for  accurate  boundary- 
layer/viscous  flow  simulation. 
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3.2  Surface  Mesh 


To  address  the  3D  tetrahedralization  problem  one  needs  to  define  surfaces  whose  union  forms 
the  boundary  of  the  domain  and  triangulate  these  surfaces.  In  this  subsection  we  describe  the 
procedures  used  to  define  the  surface  and  the  surface  triangulation  method  implemented  in 
TETMESH. 

3.2.1  Geometric  Modeling 

A  grid  generator  is  essentially  useless  without  a  powerful  technique  for  defining  the  geometry  to 
be  discretized.  TETMESH  allows  for  the  surface  geometry  to  be  specified  in  either  trimmed 
surface  or  network  format.  The  network  format  is  a  collection  of  quadrilateral  topology  of 
wireframes  with  rij  x  mi  points  defining  their  shape.  Networks  are  joined  with  other  networks 
along  common  edges  to  define  the  domain  boundary.  A  network  may  have  one  or  more 
networks  joined  to  each  of  the  edges,  but  no  edge-to-edge  overlap  is  allowed  as  shown  in 
Figure  3.1.  It  should  also  be  noted  that  the  network  may  have  degenerate  edges  (i.e.,  one  of 
the  edges  can  be  collapsed  to  a  point  to  form  a  triangular  surface  Figure  3.2). 

Trimmed  surfaces  are  defined  as  two  pieces  of  information:  a  general  surface  shape  de¬ 
scription  (the  defining  support  surface),  and  a  set  of  one  or  more  boundary  trimming  curves 
which  can  be  broken  into  multiple  edges  with  common  end  points  (Figure  3.3).  The  trimming 
curves  define  the  actual  surface  shape  and  allows  for  arbitrary  topology  surfaces  (holes,  non¬ 
convexity  etc.,).  TETMESH  converts  the  network  format  files  to  the  trimmed  surface  format 
automatically. 

3.2.2  Surface  Triangulation 

Surfaces  to  be  triangulated  are  defined  as  trimmed  surfaces  using  either  wireframe  support 
surfaces  and  point-defined  trimming  curves,  or  using  IGES  entities  126/128  and  141/143. 
These  trimmed  surfaces  are  discretized  in  two  steps:  a)  first  the  user  defines  the  spacing  on  each 
edge  (trimming  curve)  of  each  surface  using  several  available  options  (uniform,  input,  cluster) 
b)  secondly,  the  user  defines  the  required  spacing  control  of  triangles  on  the  surface  (flatness 
measure,  curvature,  spacing),  c)  finally,  the  user  requests  surface  triangulation.  Surfaces  can 
be  re-triangulated  if  the  spacing  or  quality  is  not  acceptable. 

Surface  triangles  are  generated  using  a  2-1/2  dimensional  Delaunay  scheme  as  follows: 

1.  Triangulate  all  the  boundary  nodes  in  parameter  space 

2.  Transform  to  physical  space 

3.  Swap  diagonals  to  ensure  that  the  Surface  Delaunay  Property  (SDP)  is  satisfied  (see 
below) 

4.  Incrementally  insert  new  nodes  into  the  triangles  with  the  largest  minimal  circumspheres 
in  turn  as  follows: 
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Common  Edge 
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Figure  3.1:  The  network  format  of  surface  patches.  In  this  figure  we  show  the  surface  definition  of  a  forebody 
an  assembly  of  network  patches. 


Collapsed  edge. 
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to  define  the  surface  boundary. 


Figure  3.3;  Trimmed  surface  generated  using  trimming  curves. 


•  The  new  node  has  (x,  y,  z)  coordinates  of  the  triangle’s  circumsphere 

•  Find  a  triangle  in  parameter  space  containing  the  new  node 

•  Create  three  new  triangles  from  the  new  node  and  the  containing  triangle 

•  Swap  diagonals  to  regain  the  SDP 

The  process  is  continued  until  the  desired  triangle  spacing  is  achieved.  This  procedure  has 
a  great  advantage  over  simply  triangulating  in  parameter  space.  Parameter  space  triangula¬ 
tions  suffer  from  the  distortions  of  mapping  to  physical  space.  Thus,  if  one  tries  to  achieve 
nice  triangulations  in  parameter  space,  there  is  no  guarantee  that  the  triangulation  will  be 
acceptable  in  physical  space.  Our  procedure  produces  Delaunay  surface  triangulations,  and 
thus  achieves  the  same  quality  benefits  as  planar  Delaunay  surface  triangulation  [16]. 

3.2.3  Sub-Dimensional  Delaunay  Property 

The  Delaunay  criterion  for  a  simplicial  mesh  is  that  no  node  lies  in  the  strict  interior  of  any 
simplex  circum-hypersphere.  This  concept  is  applicable  to  any  space  dimension.  The  circum- 
hypersphere  is  defined  as  the  hypersphere  in  i?”  that  passes  (uniquely)  through  the  xi  -1-  1 
vertices  of  the  n-simplex.  The  hypersphere  is  unique  up  to  degeneracies  that  occur  due  to 
pathological  placement  of  nodes  (e.g.,  all  nodes  in  an  m-dimensional  hyperplane,  vti  <  n). 
One  can  define  a  locus  of  sub-dimensional  hyperspheres  for  any  dimension  m  <  n,  which  we 
term  the  set  of  m-sphere’s  for  a  given  n-simplex.  For  example,  a  line  segment  (1-simplex) 
embeded  in  2D  has  a  set  of  1-spheres  that  have  their  centers  along  a  line  that  passes  through 
the  mid-point  of  the  line  segment  (the  center  of  the  1-sphere  in  the  line-segment  s  coordinate 
space)  and  perpendicular  to  the  segment. 

Theorem  1:  An  m-simplex  will  appear  in  an  n-dimensional  Delaunay  triangulation  if  and 
only  if  there  exists  an  m-sphere  of  the  m-simplex  that  does  not  contain  any  vertices  of  the 
triangulation. 

We  can  now  state  the: 

Surface  Delaunay  Property  (SDP):  a  surface  triangulation  is  Delaunay  if  its  trian¬ 
gles’  minimal  circumspheres  do  not  contain  other  nodes. 

The  minimal  circumsphere  of  a  triangle  embedded  in  3D  is  just  the  sphere  passing  through 
the  three  nodes  with  center  at  the  center  of  the  circumcircle  of  the  triangle  lying  in  the  plane 
of  the  triangle  (alternatively,  we  see  that  it  is  the  minimal  radius  sphere  passing  through  the 
three  triangle  vertices). 

The  SDP  gives  a  concrete  criterion  for  determining  which  surface  triangulations  are  com- 
patible  with  a  3D  tetrahedral  Delaunay  grid  computed  from  the  surface  nodes.  It  is  not 
surprising  that  most  of  the  difficulty  that  many  researchers  have  had  with  3D  constrained 
Delaunay  meshing  is  due  to  the  fact  that  most  surface  triangulations  are  simply  incompatible 
with  the  volume  mesh. 
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Figure  3.4:  Edge  swapping  criterion  for  a  pair  of  triangles,  a)  deos  not  satisfy  the 
Delaunay  criterion  b)  satisfies  Delaunay  criterion. 


Note  that  we  use  the  minimal  circumsphere  as  the  criterion  for  determining  the  Delaunay 
property  of  a  surface  triangulation.  This  is  only  a  sufficient  condition,  however,  and  simply 
states  that  if  the  minimal  condition  is  satisfied,  the  mesh  will  be  Delaunay,  and  the  surface 
triangles  will  appear  in  the  volume  mesh. 

3.2.4  Diagonal  Swapping 

We  achieve  a  Delaunay  surface  mesh  (from  a  non-Delaunay  one)  by  swapping  diagonals  on 
the  surface  for  each  pair  of  triangles  sharing  an  edge  that  contain  the  opposite  node  in  the 
triangle’s  circumsphere  center.  Edge  swapping  consists  of  modifying  an  existing  triangulation 
by  re-connecting  nodes.  In  two-dimensions,  pairs  of  triangles  are  examined  and  the  resulting 
quadrilateral’s  diagonals  are  swapped  if  the  result  produces  a  Dulanay  mesh  (see  Figure  3.4). 

3.2.5  Extrusion 

Once  the  surface  triangles  have  been  generated,  the  surface  can  be  extruded  by  user  specified 
criteria  to  produce  a  layer  of  prismatic  elements  on  the  surface  boundary.  This  layer  consists 
of  ni  sub-layers  of  prisms  that  produce  a  semi-structured  mesh  suitable  for  accurate  viscous 
flow  computations  in  the  near-body  region.  Clustering  functions  can  be  applied  to  achieve  a 
desired  first-layer  spacing  off  the  body. 

Nodes  on  the  extruded  components  are  generated  in  a  marching  scheme  away  from  the 
body.  The  scheme  is  given  as  follows: 
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1)  Compute  the  surface  normals  at  each  node. 

2)  Smooth  the  surface  normals  by  averaging  with  each  neighboring  node  (from  the  connec¬ 
tivity  graph  of  the  surface  triangulation) 

—  smoothing  is  user  controlled  by  the  number  of  passes  and  a  relaxation  factor 

3)  Create  each  layer  away  from  the  body  from  the  previous  one: 

a)  new  point  location  is  given  by  Xj+i  =Xi  +  di*  rii  where  i  is  the  layer  number,  di  is 
the  layer  spacing  and  Ui  is  the  normal  for  layer  i  at  the  given  node 

b)  After  computing  the  new  node  location,  the  new  normal  is  computed  by  smoothing 
it  from  its  neighbors  as  was  done  for  the  initial  surface  triangulation  normals. 

4)  after  nodes  have  been  created,  create  the  prisms  in  the  layer 

This  scheme  is  augmented  to  allow  for  the  addition  of  new  nodes  due  to  face  enforcement 
or  adaptivity  (see  below).  Figure  3.5  shows  prismatic  elements  extruded  from  a  triangulated 

surface  of  a  sphere. 


3.2.6  Zero-thickness  Surfaces 

When  modeling  fluid  flow  problems,  we  commonly  encounter  thin  mounted  surfaces  such  as 
fins  on  a  missile.  In  certain  instances  the  top  and  bottom  surfaces  of  such  fins  may  be  so 
close  to  each  other  that  it  would  be  beneficial  (without  loss  of  flow  solution  accuracy)  from 
the  user  perspective  to  represent  these  as  a  single  surface.  In  TETMESH  we  call  such  surfaces 
a  zero-thickness  surfaces  and  represent  these  as  a  single  trimmed  surface.  For  each  zero¬ 
thickness  surface,  after  triangulation,  TETMESH  creates  an  new  surface  with  identical  surface 
nodes  and  triangle  distributions  such  that  the  normals  are  pointed  in  opposite  direction.  These 
newly  created  surfaces  are  called  buddy  surfaces.  These  two  collapsed  surfaces  form  the  top 
and  bottom  boundaries  of  the  zero-thickness  surface  and  prevent  any  fluid  interaction  across 
the  boundary.  An  example  of  the  zero-thickness  surface  in  shown  in  Figure  3.6a.  In  this 
example  the  flap  (called  ThinJFin)  coming  out  of  the  cube  is  modeled  as  a  zero-thickness 
surface.  Figure  3.6b  shows  the  surface  normals  for  the  original  surfaces  and  the  newly  created 
buddy  surface.  Figure  3.6c  shows  the  surface  triangulation  and  the  surface  normals  for  the 
nodes  on  the  original  surface  and  the  new  surface  created  by  TETMESH. 


3.3  Volume  Mesh 

The  volume  mesh  is  computed  from  the  positions  of  the  (possibly  extruded)  surface  nodes 
and  triangles.  We  first  tetrahedralize  the  space  between  the  inner  and  outer  surfaces  using 
the  scheme  described  below.  This  forms  a  tetrahedralization  of  the  boundary  nodes  which 
then  can  be  checked  to  see  if  any  surface  triangles  need  to  be  enforced.  The  enforcement  of 
any  missing  triangles  is  performed  as  discussed  below.  Once  all  surface  triangles  exist  in  the 
tetrahedralization,  interior  nodes  are  added  at  tetrahedra  circumsphere  centers  as  was  done 
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Original  wireframe  representation  of  the  surface 
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Figure  3.5;  Illustration  of  prismatic  elements  created  by  extruding  the  surface  triangles  normal  to  the  surface. 


for  the  surface  triangulation.  We  add  nodes  in  the  tets  that  have  the  largest  circumsphere 
radius  until  the  desired  spacing  is  achieved.  Currently,  tetrahedral  spacing  is  controlled  by  a 
linear  interpolation  of  the  surface  triangle  spacing  into  the  interior  (defined  by  the  boundary 
tetrahedralization),  although  it  would  be  quite  easy  to  incorporate  a  user-supplied  spacing 
function.  However,  we  have  found  that  with  adaptivity,  it  is  only  necessary  to  achieve  a  rea¬ 
sonably  good  mesh  to  start  with,  and  the  automatic  method  for  spacing  control  in  TETMESH 
suflBces  for  that  purpose. 

3.3.1  Volume  Tetrahedralization 

We  use  a  modified  version  of  point  insertion  followed  by  diagonal  swapping  extended  to  3D. 
The  method  is  modified  from  the  original  as  described  in  [2,  15].  The  basic  idea  is  to  find  the 
tetrahedron  that  contains  the  new  node  in  its  strict  interior,  form  four  new  tetrahedra  in  the 
interior  of  that  tet  (by  joining  the  new  node  to  the  four  faces  of  the  tet),  then  performing  three- 
dimensional  edge  swapping  (3DES).  3DES  is  an  implementation  of  Lawson’s  observation  [17] 
that  there  are  only  2  unique  triangulations  of  any  configuration  of  n  -f  2  points  in  n— D  space. 
Thus,  in  3D,  two  tets  sharing  a  face  consist  of  five  vertices,  and  the  opposite  vertices  can 
be  connected  to  form  an  alternate  tetrahedralization  of  the  five  nodes.  Moreover,  three  tets 
sharing  an  edge  can  (sometimes)  be  swapped  back  to  form  two  tets  sharing  a  common  face. 
Lawson  also  showed  that  there  is  only  one  unique  triangulation  that  satisfies  the  Delaunay 
property  (modulo  degeneracies).  The  implementation  of  this  edge  swapping  algorithm  requires 
the  following: 

1)  2  — V  3  swap:  swap  two  tets  sharing  a  face  into  three  tets  sharing  a  common  edge;  this 
can  only  be  performed  if  the  two  tets  form  a  convex  polyhedron 

2)  3  — >  2  swap:  swap  three  tets  sharing  an  edge  to  form  two  tets  sharing  a  face;  this  can 
only  be  performed  if  the  resulting  polyhedron  is  convex 

3)  4  — >  4  swap:  four  tets  sharing  a  common  edge  are  swapped  to  form  four  new  tets;  this 
case  must  be  implemented  to  deal  with  the  case  where  the  polyhedron  of  cases  1)  and 
2)  is  degenerately  convex. 

4)  2  — )■  2  swap:  two  tets  on  a  boundary  sharing  a  common  edge  are  swapped  to  give  two 
new  tets. 

3.3.2  Surface  Classification 

We  can  identify  which  triangles  will  not  appear  in  the  Delaunay  mesh  using  Theorem  1  above 
as  follows.  First,  we  classify  all  triangles  whose  minimal  circumspheres  don’t  contain  other 
nodes  as  type  A.  These  triangles  will  appear  automatically  as  faces  of  the  volume  mesh.  The 
rest  of  the  triangles  are  of  type  B.  Now,  we  further  classify  the  type-B  triangles  as  follows: 

•  If  there  exists  a  2-sphere  for  the  triangle  that  does  not  contain  another  node,  classify  it 
as  B.l,  else  classify  it  as  B. 2 
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•  Furthermore,  for  B.2  triangles,  we  classify  them  as: 

B.2.1:  the  forming  points  of  the  triangle  are  not  violators  of  any  other  surface  triangle. 

B.2.2:  one  or  more  forming  point  of  the  triangle  is  a  violator  of  another  surface  triangle. 

Using  this  classification,  we  insert  each  node  of  the  surface  triangle-by-triangle  until  they 
are  all  inserted  as  follows: 

1)  Insert  all  type-B.2.2  triangles;  insert  triangle  edges  and  faces  into  the  face  checker’s  hash 
table. 

2)  Insert  all  other  type  B  triangles;  insert  edges  and  faces  into  the  face  checker  hash  table. 

3)  Insert  all  non-inserted  nodes  with  the  face  checker  ON. 

The  face  checker  uses  a  map  between  triangle  edges  (pairs  of  vertices)  and  tet  faces  (triples 
of  vertices)  to  ensure  that  once  an  edge  or  face  of  a  required  surface  triangle  is  in  the  tetrahe- 
dralization,  it  will  not  be  removed  by  any  subsequent  addition. 

Note  that  he  B.2.2  surface  triangles  are  the  most  difficult  to  enforce  since  their  nodes 
potentially  violate  the  Delaunay  property  of  other  surface  triangles. 

The  tetrahedralization  process  can  be  modeled  by  simulating  the  addition  of  the  nodes 
in  the  order  given  above.  This  simulation  then  gives  the  user  an  estimate  before  the  tetrahe¬ 
dralization  of  how  many  unenforced  faces  will  exist.  Additional  research  is  needed  on  further 
utilizing  the  triangle  classification  for  finding  an  optimal  insertion  order  within  a  particular 
sub-type.  We  speculate  that  an  ordering  can  be  found  for  triangle  insertion  that  will  further 
minimize  the  number  of  unenforced  faces  that  have  to  be  enforced. 

3.3.3  Face  Enforcement 

The  above  procedure  for  triangle  classification  and  ordered  insertion  can  still  lead  to  situations 
where  the  volume  mesh  does  not  contain  some  of  the  required  surface  triangles.  We  enforce 
these  remaining  triangles  in  a  two  step  process. 

As  noted  by  Hazlewood  [9,  10,  11]  (and  later  by  Weatherill  [24]),  one  can  enforce  the 
triangles  by  finding  the  intersection  of  the  triangles  edges  and  face  with  the  existing  mesh, 
then  inserting  new  nodes  at  the  locations  of  the  intersections.  Here’s  the  basic  algorithm: 

1)  Enforce  all  three  edges  of  the  required  surface  triangle.  If  triangle  edge  E  does  not 
already  exist  (possibly  with  additional  vertices  along  it)  then  that  edge  pierces  at  least  one 
tetrahedron.  Subdivide  each  pierced  tetrahedron  to  create  the  segment  of  edge  E  passing 
through  it. 

2)  Enforce  the  required  triangle  face.  If  face  ABC  does  not  already  exist  (possibly  contain¬ 
ing  vertices  on  its  boundary  that  were  created  in  step  1)  then  that  face  intersects  at  least  one 
tetrahedron.  Subdivide  each  intersected  tetrahedron  to  create  the  portion  of  face  ABC  that 
intersected  it.  Note  that  the  portion  of  ABC  that  is  common  to  an  intersecting  tetrahedron 
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may  be  triangular  or  quadrilateral.  The  union  of  these  portions  of  triangle  ABC  then  form 
the  required  face  to  be  enforced. 

In  the  above,  each  intersected  tetrahedron  is  subdivided  into  smaller  tetrahedra  with  the 
constraint  that  the  specified  edge  or  face  be  a  part  of  the  resulting  sub-tetrahedralization. 
When  there  is  more  than  one  possible  sub-tetrahedralization,  a  Delaunay  one  is  chosen. 

The  next  step  is  to  remove  newly  added  nodes  by  collapsing  edges  between  such  nodes  and 
other  nodes  in  the  enforced  face. 

1)  for  each  new  face  (interior)  node,  try  to  collapse  it  to  another  face  node,  another  edge 
interior  node,  or  a  node  of  the  original  face. 

2)  for  each  edge  interior  node,  collapse  the  node  to  one  of  its  neighbors  on  the  edge. 

The  collapsing  procedure  is  described  below.  The  reason  that  we  try  to  remove  the  new 
nodes  added  by  the  face  enforcer  is  that  they  can  be  arbitrarily  close  to  other  nodes,  and  thus 
produce  arbitrarily  small  (and  thus  bad)  triangles  and  tets. 

3.3.4  Edge  Collapsing  Procedure 

A  major  contribution  of  this  project  has  been  the  development  of  the  following  algorithm 
for  collapsing  edges  to  remove  unwanted  nodes.  Other  schemes  for  node  removal  require  re- 
tetrahedralizing  the  cavity  that  results  from  node  removal  [24].  It  is  well  known  that  arbitrary 
cavities  cannot  be  tetrahedralized  while  maintaining  the  boundaries.  This  is  just  another 
example  of  the  difficulty  of  3D  tetrahedral  mesh  generation.  Moreover,  merely  determining 
whether  a  configuration  is  tetrahedralizable  has  been  shown  to  be  NP-Complete  [19].  However, 
the  following  constructive  algorithm  provides  both  a  test  for  tetrahedralizability  as  well  as  the 
resulting  tetrahedralization,  in  the  case  that  one  starts  with  an  existing  tetrahedralization,  and 
removes  a  given  node. 

1)  for  each  possible  collapse,  check  if  the  resulting  tetrahedralization  is  valid:  each  tet  that 
doesn’t  have  a  collapsed  edge  must  have  positive  volume  2)  continue  until  a  valid  collapsed 
configuration  is  found  3)  remove  all  tets  that  have  a  collapsed  edge,  and  also  the  newly  inserted 
node. 

This  edge-collapsing  procedure  is  also  used  in  STORESIM’s  adaptive  unrefinement  strategy 
(see  below).  The  above  procedure  could  be  modified  to  not  only  choose  a  valid  configuration, 
but  to  choose  the  best  quality  tetrahedralization. 

3.3.5  Boundary  Tet  Layer 

As  is  done  in  the  Rampant  code  [3],  we  introduce  a  layer  of  nodes  right  at  the  boundary 
before  face  enforcement  to  produce  an  optimal  set  of  tets  near  the  boundary.  The  nodes  are 
placed  along  a  normal  to  each  face  circumcenter  at  a  distance  of  r-\/2  where  r  is  the  minimal 
circumsphere  radius.  This  initial  layer  produces  good  quality  tets  right  at  the  boundary  and 
furthermore  reduces  the  possibility  of  badly  shaped  triangles/tets  being  generated  by  the  face 
enforcer.  No  edge  or  face  of  an  existing  surface  triangle  is  removed  by  this  process  since  the 
’face  checker’  is  on  during  the  node  addition. 
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3.3.6  Mesh  Quality 


We  define  the  aspect  ratio  of  a  tetrahedron  by  the  radius  of  its  inscribing  sphere  divided  by  the 
radius  of  the  circumscribing  sphere.  This  measure  is  scaled  so  that  an  equilateral  tetrahedron 
has  aspect  ratio  1,  and  a  zero- volume,  degenerate  tet  -  or  sliver  -  has  aspect  ratio  0. 

Slivers  are  removed  by  a  process  in  which  an  edge  of  the  degenerate  tet  is  deleted.  This 
edge  removal  is  precisely  the  3 — >2  swap  in  its  simplest  form,  but  in  general  it  is  an  n — >-2n— 4 
swap.  It  is  described  as  follows.  For  each  edge  of  the  tet  we’re  trying  to  remove,  consider 
all  tetrahedra  that  share  that  edge.  Suppose  the  candidate  edge  is  surrounded  by  n  tets. 
Those  tets  are  defined  by  n  H-  2  nodes,  two  of  which  define  the  shared  edge.  Enumerate  all 
possible  triangulations  of  the  n  nodes  not  on  the  shared  edge.  The  triangulation  takes  into 
account  no  geometry  other  than  the  cyclic  ordering  of  the  nodes.  That  is,  for  the  purpose  of 
enumerating  the  triangulations,  one  could  assume  that  the  n  nodes  lay  on  a  circle.  For  each  of 
those  triangulations  (consisting  of  p  =  n  -  2  triangles),  we  form  the  resulting  2p  tetrahedra  by 
joining  each  triangle  with  the  two  nodes  of  the  shared  edge.  This  obviously  removes  the  edge. 
Of  the  possible  resulting  tetrahedralizations,  we  choose  the  one  that  maximizes  the  minimum 
aspect  ratio  of  the  2p  new  tetrahedra.  It  should  be  mentioned  that  although  this  procedure 
is  usually  effective,  it’s  not  always  possible  to  remove  an  edge  at  all,  and  even  if  it  is  possible, 
doing  so  may  produce  worse  tets.  Also,  due  to  the  exponential  complexity  of  enumerating  all 
triangulations,  this  procedure  is  effective  only  for  very  small  n,  say  n  <  10. 

3.4  Adaptivity 

For  accurate  simulation  of  both  steady  and  unsteady  flows,  one  must  use  an  adaptive  meshing 
to  capture  important  flow  features  such  as  shock  waves  and  boundary  layers.  A  critical 
component  of  an  adaptive  code  is  the  ability  both  to  refine  and  to  unrefine  the  mesh  as 
outlined  below. 

3.4.1  Refinement 

Adaptive  refinement  of  the  mesh  is  accomplished  by  adding  new  nodes  on  selected  edges  of 
the  mesh.  This  creates  2n  new  tets  where  n  is  the  number  of  tets  surrounding  the  broken 
edge.  The  edges  are  selected  for  refinement  based  on  typical  gradient  error  indicators  (such 
as  the  density  gradient). 

3.4.2  Unrefinement 

Unrefinement  is  implemented  in  STORESIM  by  identifying  nodes  to  be  removed.  This  is 
accomplished  by  marking  all  nodes  as  NOT-REMOVED  then  marking  which  edges  have  error 
indicators  smaller  than  the  UWREFINE-CUTOFF  value.  Each  edge  then  marks  its  two  nodes,  and 
if  all  edges  surrounding  a  node  should  be  unrefined,  then  the  node  is  marked  for  unrefinement. 
Next,  the  node  is  collapsed  to  one  of  its  neighbors  by  selecting  the  neighbor  that  produces  the 
best  resulting  mesh.  Note  that  this  procedure  requires  no  cavity-filling  as  do  other  schemes  [8]. 


18 


Also,  no  expensive  tree  structure  of  tets  needs  to  be  implemented,  since  whenever  the  mesh 
is  refined,  the  old  tets  are  simply  discarded.  The  old  tets  are  easily  recovered  by  the  edge 
collapsing  procedure  (unrefinement)  as  described  above. 

3.4.3  Surfaces 

When  a  new  node  is  added  on  an  edge  during  refinement,  it  is  possible  that  the  edge  lies 
on  the  surface  boundary.  Thus,  we  must  update  the  surface  triangulation  (creating  four  new 
surface  triangles  from  the  two  sharing  the  edge)  as  well  as  ensure  that  the  new  node  lies  on 
the  actual  surface  and  not  simply  on  the  straight  line  between  the  two  edge  endpoints.  This 
ensures  high  surface  fidelity. 

Surfaces  that  have  extruded  prisms  must  have  the  prisms  adapted  as  well.  We  produce  4 
new  prisms  for  each  original  2  in  each  layer  of  the  extrusion,  and  ensure  that  the  node  on  the 
interior  boundary  lies  on  the  actual  surface. 

3.5  Interface  to  Other  Codes 

Once  the  mesh  generation  is  completed,  TETMESH  allows  the  user  to  save  the  mesh  in  several 
different  forms.  This  enables  the  user  to  output  the  mesh  and  use  it  with  a  third-party  flow 
solver.  At  the  present  time,  TETMESH  is  capable  of  outputting  the  mesh  in  the  following 
formats: 

i)  Internal  binary /ascii  format  to  be  read  into  TETMESH, 

ii)  Gamma  format:  to  input  data  to  stand-alone  STOREFLOW  solver, 

iii)  FAST  format, 

iii)  COBALT  format,  and 

iv)  ascii  format  listing  tets  and  their  faces. 

Of  these,  the  FAST  format  is  widely  published. 
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4 


Flow  Solver 


4.1  Background 

STOREFLOW  was  developed  within  the  context  of  the  STORESIM  project,  whose  primary 
goal  was  the  numerical  computation  of  the  trajectories  of  stores  (i.e.,  missiles,  bombs,  etc.,) 
released  from  flying  aircraft.  Any  solver  employed  for  such  a  purpose  must  meet  several 
stringent  requirements.  Firstly,  it  should  be  capable  of  modeling  very  complex  and  arbitrary 
geometries  that  characterize  “stores”  and  aircraft.  It  must  compute  physically  relevant  flow 
solutions  over  these  configurations  over  a  wide  range  of  speed  regimes,  from  low  subsonic  to 
transonic,  supersonic  and  hypersonic  Mach  numbers.  Finally,  it  must  deliver  solutions  for 
multiple  bodies,  moving  at  varying  relative  speeds.  This  last  capability  differentiates  it  from 
traditional  flow  solvers  that  are  mainly  used  to  compute  flows  at  some  fixed  Mach  number. 

Keeping  the  above  requirements  in  mind,  the  algorithm  used  to  develop  STOREFLOW 
was  based  on  the  SUPG  (Streamline  Upwind  Petrov  Galerkin)  method  [14]  ,  [22],  which  has 
in  recent  years  emerged  as  a  very  robust  approach  for  Finite  Element  Computational  Fluid 
Dynamics.  STOREFLOW  contains  several  new  features  to  enhance  its  capabilities.  Among 
the  most  important  in  this  list  are: 

•  A  modification  [23]  to  allow  ALE  (Arbitrary  Lagrangian  Eulerian  )  computations.  This 
essentially  makes  it  possible  to  use  moving  meshes  to  model  multiple  moving  bodies. 
This  new  ALE  formulation  is  extremely  robust,  and  can  even  produce  accurate  results 
for  bodies  launched  at  high  speed  into  ambient  air. 

•  A  very  sophisticated  far  stream  boundary  condition  designed  to  prevent  spurious  energy 
build  up  and  consequent  reflections  at  transonic  speeds.  This  approach  also  includes 
modifications  to  account  for  moving  mesh  boundaries. 

•  Euler  and  Navier-Stokes  versions.  Currently,  the  viscous,  Navier-Stokes  code  in  oper¬ 
ational  only  for  laminar  flow.  Turbulence  modeling  based  on  a  simple  zero  equation 
model  has  been  tested  for  two-dimensional  flows.  The  extension  to  three-dimensional 
turbulence  modeling,  although  straightforward,  is  not  yet  operational. 

•  Real  gas  capability  [5],  based  on  a  5  species  (Nitrogen,  Oxygen,  Atomic  Nitrogen,  Atomic 
Oxygen  and  Nitric  Oxide)  model  for  dissociated  air. 

•  Fully  integrated  with  a  6  degree  of  freedom  (6  dof)  solver  for  rigid  body  dynamics,  and 
a  mesh-mover,  that  moves  the  mesh  in  response  to  body  motions. 

•  Multiple  element  types  (tetrahedra,  Hexahedra  or  Prisms).  Additional  element  types 
can  be  easily  incorporated  into  the  code. 

4.0 
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4.2  Governing  Equations 

The  equations  governing  the  flow  of  a  compressible,  viscous  fluid  are  the  Navier-Stokes  equa¬ 
tions  and  are  given  by  the  vector  relation: 


U,,  +  =  Ff, 


where 


U  :  vector  of  conservation  =  p{ 


1 

Ui 

U2 

Us 

e 


p:  fluid  density 

Uji  j  =1,2,3  =  fluid  velocity  components 
e:  total  energy  per  unit  mass  =  Cv9  -h  ^UiUi 
9:  temperature  of  the  fluid 
C^:  specific  heat  at  constant  volume 
Fi'.  Euler  flux  vector  and  is  given  by: 


Fi 


] 


UiU  -1-  p  < 


^21  > 
^3i 


Ff  :  diffusive  flux  vector  and  is  given  by: 


Ff=< 


(  0 

Tli 
T2i 
Tsi 


> 


^  TijUj  Qi  J 

Tij  :  is  the  stress  tensor  =  Xuk,kSij  -I-  p  {uij  +  Uj^i) 

p  :  viscosity  coefficient 

A:  bulk  viscosity  coefficient  —  ~p 

Qi  :  heat  flux  component  =—k9^i 

p:  pressure 


(4.1) 


(4.2) 


(4.3) 
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4.2.1  Euler  Equations 


Consider  an  inviscid,  adiabatic  flow  field.  The  equations  governing  the  flow  of  such  fluid  can 
be  obtained  from  the  Navier-Stokes  equations  by  neglecting  all  shear  stress  terms  and  heat 
conduction  terms.  The  resulting  set  of  equations  are  called  Euler  equations.  Thus,  the  Euler 
equations  are: 


U,t  +  Fi,i  =  0  (4.4) 

The  Euler  equations  are  used  for  flows  at  high  Reynolds  numbers.  From  Prandtl’s  boundary 
layer  analysis,  this  is  a  valid  approximation  for  such  flows  outside  the  viscous  regions  close  to 
the  surfaces. 

4.2.2  Navier  -  Stokes  Equations  in  an  ALE  (Arbitrary  Lagrangian  Eulerian)  Ref¬ 
erence  Frame 

In  this  section  we  present  the  Navier-Stokes  equations  in  an  ALE  reference  frame.  The  Navier- 
Stokes  equations  as  presented  in  (4.1)  are  only  valid  in  a  fixed  frame  of  reference.  This  is  the 
conventional  “Eulerian”  framework.  Thus,  the  flow  solution  is  defined  at  arbitrary  points  “p”, 
that  are  stationary  with  respect  to  an  inertial  frame  of  reference  [x,y,z]  (see  Figure  4.1). 

In  the  “Lagrangian”  representation,  the  point  “p”  has  a  velocity  uf,  that  is  exactly  the 
same  as  the  fluid  velocity  lij.  However,  in  an  ALE  approach,  uf  is  completely  arbitrary.  It  is 
advantageous  to  use  such  an  ALE  representation  to  model  the  flow  over  multiple  bodies  that 
are  moving  at  arbitrarily  different  speeds.  The  Navier-Stokes  Equations  for  the  ALE  frame 
can  be  shown  to  be: 


(4.5) 


The  first  term  on  the  R.H.S  of  (4.5)  is  a  grid  dilation  term.  The  flux  vector  Fi,  is  defined 
as: 


Fi  =  F^-  ufU  (4.6) 

where  F^  has  the  same  form  as  the  flux  vector  in  an  Eulerian  frame.  It  has  been  shown 
[23]  ,  that  the  grid  dilation  term  may  be  conveniently  eliminated  to  give  the  ALE  equation  in 
the  form: 


where 


U,t  +  [Ai-ufI]U,i 


(4.7) 
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Ai  —  =  Euler  Flux  Jacobian  (4.8) 

Equation  (4.7)  forms  the  basis  of  the  STOREFLOW  code. 

4.3  Boundary  Conditions 

Symmetry  and  specified  heat  fiux  conditions  are  imposed  as  natural  boundary  conditions,  by 
merely  computing  the  associated  boundary  integrals.  The  “no-slip”  condition  on  solid  walls 
is  imposed  as  the  “Dirichlet”  conditions,  Ui  =  uf.  At  farstream  boundaries  the  formulation 
applies  specified  farstream  conditions  in  conjunction  with  a  non-reflective,  characteristic  theory 
[21]. 

4.4  Weak  Formulation  of  the  ALE  Navier-Stokes  and  Euler  Equa¬ 
tions,  in  the  Context  of  the  SUPG  Method 

The  weak  form  of  these  equations  ,  in  the  starting  point  for  a  Finite  Element  discretization 
and  solution  algorithm.  The  intent  here  is  to  present  the  weak  form  in  sufficient  detail  to 
facilitate  an  in  depth  examination  of  the  corresponding  computer  code.  Precise  details  on 
steps  required  to  arrive  at  this  formulation  may  be  found  in  Reference  [14], [22],  and  [23]. 
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For  the  ALE  Navier-Stokes  equations  we  employ  the  weak  form: 


/  W-{U,t  +  [Ai-ulI]ui)dn+  f  W,i-Ff-dn-  f  W-iFi-Fioo)nidr 

Jn  Jn  Jtoo 

—  j  W '  {Fi  —  Fip)  Tiidr  —  j  W  •  Fidni  •  dT  +  Isupg  +  Ishock  —  0 

•'Twall  •'^Wall 

where 

'  0  ' 

Sii 

Fip  =  P{  S2i  ^  (4.9) 

53: 

•  W  =  vector  of  weighting  functions  chosen  from  a  space  of  compact  support 

•  Isupg  =  Streamline  Operator 

•  Ishock  =  Shock  Operator 

The  streamline  and  shock  operators  are  defined  in  Appendix- A. 

4.5  Finite  Element  Discretization 

A  Finite  Element  discretization  based  on  linear  elements  has  been  implemented  for  various 
element  types.  The  elements  considered  are  triangle  and  quadrilateral  in  two-dimensions,  and 
tetrahedra  and  prisms  in  three-dimensions.  A  semi-discrete  formulation  leads  to  the  following 
representation  for  the  time  derivative  of  the  solution  variable  at  each  node,  at  time  t 

[M"]  =  -FT  (4.10) 

where  M”  is  a  mass  matrix  assembled  as  follows: 

M"  =  AN4[M"]  (4.11) 

where  Nel  —  number  of  elements, 

M"=  ,p,g  =  l,2,...,Nen  (4.12) 
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where  N^n.  =  number  of  elements  nodes, 


M‘,  =  /  iVXW<iS2  (4.13) 

where  I  is  the  5x5  unit  matrix  and  AT®,  JV®  are  linear  element  shape  functions.  Also, 
u  =  5  X  Nnp  vector  of  conservation  variable  at  all  nodes. 

Finally,  iJ”  is  a  5  x  Nnp  residual  vector,  which  we  assemble  as  follows: 


(4.14) 

R®=  [R;],p  =  l,2,...,A-en 

(4.15) 

II 

n;  .  [AoVt  -h  if dD  +  {Isupg  +  hkock)  dn 

infty 

Fioo)  riidT  -  f  Nf  {Fi  -  Fip)  UidT  -  f  N^  '  Ff  ■  UidT 

Jt‘  ^  Jr' 

(4.16) 

where 


V  =  vector  of  entropy  variables 
Ao,  =  flux  matrices 


) 


Defined  in  the  Appendix  —  A 


4.5.1  Time  Discretization 

The  semi-discrete  form  (4.10)  has  been  used  to  developed  an  explicit  Runge-Kutta  and  implicit 
schemes  of  varying  orders  of  accuracy.  For  example,  the  general  “a”  scheme  for  computing 
the  solution  at  the  (n  -I-  1)**  time  step  is; 


+  At  -1-Fa  u”  j  0  <  a  <  1 


(4.17) 


if 

a  =  0  :  Explicit^  Forward  Euler 
Oi  =  1  :  Implicit,  Backward  Euler 
a  =  i  ;  2^^ order  accurate,  Cranck  —  Nicholson  scheme 
In  a  similar  manner,  Runge-Kutta  schemes  of  varying  orders  (3  ~  5)  have  been  coded. 
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4.5.2  Non-Linear  GMRES,  IMPLEX-Solver 


An  Implex  (Implicit /Explicit)  solver,  for  the  Cranck-Nicholson  scheme  has  been  developed. 
This  automatically  picks  Implicit  and  Explicit  nodes  in  the  mesh  based  upon  a  user  supplied 
range  for  the  CFL  coordinates  as  follows: 

•  Choose  CFLmin  suid  CFL  max 

•  Compute  Ate  =  element  “e”  ,  where 

he  =  element  length  parameter 
V  =  velocity 
a  =  Speed  of  Sound 

•  Compute  Ati  for  each  node  as  the  minimum  value  of  the  elements  connected  to  that 

node.  Then  find  Atmin  =  min(Ati, . ,  AU),  which  is  the  timestep  used  to  advance  the 

solution. 

•  For  each  node  compute  CFLi  =  cflZZ  ' 

•  If  CFLi  is  less  than  CFL^n 

Node  “i”  is  explicit 
else 

Node  “i”  is  implicit 
endif 

Both  the  fully  implicit  Backward  Euler  and  the  implicit  algorithm  require  the  solution  of 
a  large  set  of  Non-linear  equations,  which  can  be  represented  as: 


N(u)  =  {0} 


(4.18) 


This  equations  has  5  •  iV„p  unknowns. 

A  Newton-Raphson  linearization  yields 


N  [vf  -h  Au^)  =  N  (n^)  + 


dN 


An*  Ri  0 


du 

—  Avf 


(4.19) 


One  needs  to  solve  this  arge  linear  system  of  equations,  for  the  Kth  approximation  to  u. 


du 


Au'^ 


(4.20) 


26 


This  is  repeated  until  Au*  becomes  very  small,  for  some  value  k.  In  practice,  (4.20)  is  very 
expensive  in  terms  of  computer  memory.  Thus,  recourse  is  available  through  the  Nonlinear 
GMRES  (Generalized  Minimum  Residual)  Method  [4]  &  [25].  In  this  method,  the  GMRES 

• 

^  •  p,  where  p  is  a 

Krylov  vector  (see  references  above).  This  product  is  approximated  to  second  order  accuracy 
in  the  Nonlinear  GMRES  procedure  as  follows. 


aivf 


du 


■P  = 


N  {u  +  ep)  —  N{u  —  ep) 

Te 


(4.21) 


4.5.3  Preconditioning 

A  preconditioner  for  (4.20)  based  on  an  approximate  inverse  to  has  been  developed.  The 

["  1  A; 

^1  .  This 

makes  the  inversion  very  fast.  Another  preconditioner  that  has  been  implemented  is  [18]. 


K  = 


(4.22) 


where  K~^  =  L  •  n  and  K  ^ 
ignored. 


is  computed  such  that  the  “fill  -in”  during  factorization  is 


4.6  Flow  Solver  Capabilities 

As  mentioned  earlier,  STORESIM  contains  a  module  STOREFLOW  that  implements  a  vec¬ 
torized  version  of  the  Streamline  Upwind  Petrov-Galerkin  (SUPG)  method  [13].  The  SUPG 
method  has  been  modified  (as  outlined  in  [23])  to  incorporate  the  extra  terms  that  arise  from 
the  mesh  movement  in  an  Arbitrary  Lagrangian-Eulerian  (ALE)  formulation  of  the  weak  form 
of  the  Navier-Stokes  equations.  Key  features  of  STOREFLOW  are: 

•  Physics: 

-  inviscid  (Euler) 

-  viscous  (Navier-Stokes) 

*  laminar 

*  turbulent  -  0-equation  model 

-  real  gas 

*  equilibrium  chemistry 
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*  air  model  with  5  species 


•  Time: 

—  steady-state 

—  time-accurate 

•  Solver: 

—  explicit 

—  implicit 

—  implicit/explicit  (IMPLEX) 

•  Mesh: 

—  static  mesh 

—  dynamic  mesh 

*  ALE  formulation 

Any  combination  of  the  above  physics,  time  stepping,  solver  or  mesh  can  be  chosen. 

STOREFLOW  has  an  element  library  that  includes  a  linear  triangle,  quadrilateral,  prism, 
tetrahedron,  hexahedron  elements.  It  has  a  unique  software  architecture  in  which  a  single 
source  code  is  maintained  which  is  pre-processed  (by  the  Unix  m4  macro  processor)  to  produce 
actual  FORTRAN  code  for  2D  or  3D,  any  element  type  (or  combination  of  types  such  as  tets 
and  prisms),  or  time-stepping  scheme  (explicit,  implicit,  IMPLEX)  as  desired. 
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5  Six-DOF  Rigid  Body  Dynamics  Solver 


STORESIM  has  a  fully  integrated  (coupled)  six  degree-of-freedom  rigid  body  dynamics  (6DOF- 
RBD)  solver  that  uses  the  consistently  computed  forces  and  moments  on  the  specified  moving 
bodies  to  accurately  predict  the  trajectory  of  the  moving  bodies  under  either  free  or  powered 
flight.  The  mass  properties  of  each  moving  body,  thrust  and  moment  forces  as  well  as  option¬ 
ally  prescribed  (captive)  trajectory  are  specified  in  a  simple  C-like  language  in  the  input  to 
STORESIM. 

The  6DOF-RBD  solver  uses  a  fully  general  quaternion  formulation  to  avoid  the  singularities 
inherent  in  the  use  of  the  Euler  angle  formulation.  This  results  in  a  system  of  14  ordinary 
differential  equations  governing  the  motion  of  the  body. 

5.1  Six-DOF  Rigid  Body  Dynamics 

5.1.1  Basic  Equations 

The  fundamental  equations  describing  the  six  degrees-of-freedom  (6  -dof)  dynamics  of  a  rigid 
body,  are  expressed  with  respect  to  an  inertial  frame  [X,Y,Z],  in  conjunction  with  a  rotating 
reference  frame  X,  Y,  Z,  that  is  fixed  to  the  body  (see  Figure  5.1).  For  store  simulations  where 
flight  durations  are  relatively  short,  an  earth  axis  system  [X,Y,Z]  is  a  good  approximation  to 
an  inertial  frame.  Conversely,  for  long  duration  motions,  it  would  be  necessary  to  include  the 
effect  of  the  earths  rotation  with  respect  to  a  star  fixed  inertial  frame.  In  the  present  study, 
it  is  assumed  that  [X,Y,Z]  is  earth  fixed. 

These  fundamental  equations,  which  are  in  fact  expressions  of  Newtons  Laws  for  the  ma¬ 
terial  particles  forming  the  body  are: 


Fext  =m*fo 

^^ext  ~  Fq 

where, 

m  =  mass  of  the  body 

ro  =  position  vector  of  the  bodys  mass  center  “0” 

Hq  =  angular  momentum  vector 

Fext  =  external  forces  arising  from  aerodynamic  and  gravitational  loads. 

Mext  =  external  moments  arising  from  aerodynamic  and  gravitational  loads. 

5.1.2  Rates  of  Change  of  Linear  and  Angular  Momentum 

Equation  (5.1)  (5.2)  relate  externally  impressed  forces  {Fext  and  Mext  )  and  inertia  forces 
developed  by  the  rigid  body.  As  these  forces  are  expressed  in  the  inertial  frame  [X,  Y,  Z], 
it  is  essential  to  evaluate  the  rates  of  change  of  linear  and  angular  momentum  as  seen  by  an 
observer  [X,Y,Z]. 


(5.1) 

(5.2) 
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Figure  5.1:  Coordinate  system  for  the  fluid  domain  and  the  moving  body. 


5.1.3  Rate  of  Change  of  Linear  Momentum 

This  is  written  as  m*fo.  It  is  convenient  to  express  the  same  in  terms  of  the  velocity  compo¬ 
nents  (So  of  the  mass  center.  The  (  ~  )  is  used  to  signify  that  So  represents  the  components 
of  the  velocity  vector  in  a  co-ordinate  system  that  is  instantaneously  aligned  with  the  body 
axes  [X,y,  Z  ].  Likewise,  Fgxt  represents  the  components  of  the  external  force  Fgxt- 

Now,  the  linear  momentum  vector  as  seen  by  an  observer  in  [X,y,  Z],  is  m  *  So-  The 
rate  of  change  of  this  vector,  as  seen  by  an  observer  in  [X,Y,Z],  is  obtained  with  the  aid  of  a 
standard  result  governing  rotating  reference  frames,  viz: 


^xyz 


—  Axyz  +  UJb  ®  A 


(5.3) 


where 

^b= 


P 

q 

r 


>  =  rotational  velocity  of  the  body. 


Thus,  using  (5.3)  we  obtain; 


muo  =  rnuoliys  -b  0)5  (gi  muo  =  F^ 


ext 


(5.4) 
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5.1.4  Rate  of  Change  of  Angular  Momentum 

This  is  determined  in  a  manner  analogous  to  that  used  for  evaluating  the  linear  momentum. 
The  angular  momentum  vector  seen  by  an  observer  in  [X,  Y,  Z]  is: 

t^o\x,y,z  —  I (^•^) 
where  I  is  the  symmetric  inertia  tensor  of  the  body  in  the  body  axes: 


Thus,  as  in  (5.3),  we  obtain 


(5.6) 


■Hoixyr  —  -^olxyz  "t"  ®  -^folxyz 

=  iu!b  +  OJb®  i^^b)  —  ^ext 

Finally  from  (5.4)  and  (5.7  )  we  get 


(5.7) 


““  1 

^ 2Xt 

-1 

Uq  =  <( 

m 

-  a>6  (8)  uo 

1 

1  ^^^0  J 

(5.8) 


=  < 


P 


i  =  ^  *  [Meit  —  Qb®  (^‘^ft)] 


(5.9) 


5.1.5  Transformation  Between  the  Body  Fixed  and  Inertial  Axes  Systems 
The  orientation  of  the  body  axes  [X,Y,  Z]  with  respect  to  [X,Y,Z]  is  needed  for  two  reasons: 

•  To  determine  the  components  of  the  external  forces  and  moments,  as  required  by  the 
6-dof  equations. 

•  To  determine  the  precise  location  of  the  desired  points  on  the  body  surface.  Note,  that 
the  6-dof  equations  calculate  only  the  position  of  the  mass  center. 
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Now,  in  principle,  the  rotational  velocity  Cbi  can  be  integrated  in  time  to  obtain  the  body’s 
orientation.  However,  it  is  necessary  to  express  any  arbitrary  rotation  in  terms  of  a  sequence 
of  3  independent  rotations.  This  is  done  in  a  standard  manner  using  “Eulerian  angles”,  that 
we  present  in  the  next  section,  where  it  will  be  seen  that  the  time  rates  of  change  of  the 
Eulerian  angles  are  in  fact  related  to  the  rotational  velocity.  Finally,  we  present  the  so  called 
quaternionic  variables  (Euler  parameters)  whose  use  is  preferred  to  the  Eulerian  angles  in 
numerical  computations  of  6-dof  dynamics. 

These  are  defined  by  the  three  elementary  rotations  ^  {yaw),  6  (pitch),  and  ^  (roll)  which 
lead  the  earth  fixed  axes  [X,Y,Z]  (see  Figure  5.2).  [Xi,  Yi,  Zi]  and  [X2,  Y2,  Z2]  are  intermediate 
positions  of  the  axes  at  each  stage  of  the  Eulerian  rotations.  These  facilitate  the  formulation 
of  the  rotation  matrix  between  the  earth  and  body  axes. 

5.1.6  Rotation  Matrices 

The  sequence  of  rotations  described  by  the  Eulerian  angles  leads  to  the  following  rotation 
matrices,  G,j„G0,G^. 


G  / 

[x,y,z]  4  [xx,yi,zi] 

[xi,yi,zi]  ^  [x2,yi,Z2] 

N,  y2,  Z2]  ^  [x,  y,  z]  (5.10) 


where 


’  cosip 

sinip 

0  ■ 

cos6 

0 

—sind 

G(p  = 

—sinip 

cosip 

0 

Ge  = 

0 

1 

0 

0 

0 

1 

sin9 

0 

cos6 

G^  = 


10  0 
0  cos(j)  sin(p 
0  —sincj)  coscj) 


Thus,  [X,  y,  Z]  4  [X,  Y,  Z]  where 


G  =  G^  ■  Gq  •  G,p 


(5.11) 


Note  that 


=  Gl  :  G;^  =  GJ 

G^^  =  GJ  :  G-^  =  G^  (5.12) 
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5.1.7  Relation  Between  Rotational  Velocity  and  the  Time  Derivatives  of  the 
Euler  Angles 

The  rotational  velocity  of  the  body,  Ub  has  components  [p,q,r]  in  the  body  axes  system.  Now, 
stemming  from  the  definition  of  the  Euler  angles,  we  can  also  express  Ub  in  terms  of  its 
components  along  the  directions  of  the  three  elementary  rotations.  Thus,  we  obtain. 


Ub  =  ‘tP  •  Tiz-^  +6y^+  ^ 


n^r 


(5.13) 


where  nzy^Uy^^ns;  are  unit  vectors  along  directions  Z1,Y2,X  respectively. 
Using  (5.13  )  in  conjunction  with  (  5.10  )  one  can  show  that 


(p  —  ipsinO  1 

'  P  ' 

ipcosB  •  sincp  +  6cos(p  , 

q  > 

,  ipcos6cos(p  —  dsin4>  J 

1  r  , 

(5.14) 


It  is  clear  from  (5.14)  that  p,q,r  can  be  integrated  with  respect  to  time  to  give 
However,  as  we  shall  see  in  the  next  section,  the  Eulerian  angles  are  not  the  most  suitable 
parameters  to  use  for  rigid  body  dynamics,  and  that  it  is  preferable  to  employ  the  so  called 
“quarternionic  variables  “(Euler  parameters)  in  their  place. 


5.1.8  Singular  Behavior  of  the  Rotational  Derivatives 

From  (5.14),  we  see  that  the  time  derivative  of  the  Euler  angles  are  related  to  the  rotational 
velocity  components  by  the  relation: 


^  P  ' 

6 

q  > 

^  J 

where 


(5.15) 


E  = 


—sin6  0  1 

cosO  •  sinp  cos4>  0 
cosO  ■  cosp  —sincf)  0 


(5.16) 


Now  Det{E)  =  —cosB  =  0  (if  0  =  ±7r/2)  Hence,  at  6  =  ±7r/2,  equation  (5.14)  is  singular. 
Geometrically,  this  stems  from  the  fact  that  at  ^  =  ±7r/2,  ip  is  undefined  and  can  take  any 
value.  The  effect  of  this  singularity  on  numerical  computations  is  to  cause  increasely  large 
errors  in  the  vicinity  oi  6  =  ±7r/2,  leading  to  completely  erroneous  results.  The  standard 
remedy  employed  in  this  case  is  to  use  “quarternionic  variables  ” . 
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5.2  Quarternions  (Euler  Parameters  ) 

Any  arbitrary  rotation  in  3-D  space  is  equivalent  to  a  pure  finite  rotation  (a)  about  a  uniquely 
defined  axes.  This  corresponds  to  the  mathematical  fact  that  the  rotation  matrix  governing  a 
rotation  has  a  single  real  eigen  value  “A”  and  a  unique  eigenvector  (with  normalized  compo¬ 
nents  rii).  This  is  equivalent  to  a  pure  rotation  of  a  =  A  along  the  axis  n  —  i. 

The  four  quarternions  are  then  defined  as  follows: 


z=  niSina/2  :  i  =  1, 2, 3 


94  =  cosa/2 


(5.17) 


It  is  easily  verified  that 


9l+92  +  ?3+9l  =  l 

The  rotation  matrix  can  be  expressed  in  terms  of  the  quarternions  as: 


(5.18) 


G  = 


■2(91  +  9|)-1 
2(9i92  -  mi) 
.  2(9i93  -i-  mi) 


2(91 92  +  mi) 

m + 9i)  - 1 

2(9293  -  9194) 


2(9193  -  9294)  ' 
2(9293  +  9194) 
2(9l  +  9|)-l  J 


(5.19) 


Also,  the  time  derivatives  of  the  quarternions  are  related  to  the  rotational  rates  by: 


< 


■  94 

-93 

92 

-9i  ■ 

92 

1 

93^ 

94 

-9i 

-92 

?3 

'  “  2 

—92 

9i 

94 

—93 

/ 

r 

94  . 

.  -9i 

-92 

-93 

-94  . 

.  0 . 

(5.20) 


Now,  it  is  evident  from  (5.20)  that  the  time  derivatives  of  the  quarternions  are  always 
uniquely  defined.  This  is  contrary  to  that  for  the  Eulerian  angles,  which  have  a  singularity  at 
6  =  ±7r/2.  Thus,  (5.20)  is  used  instead  of  (5.14)  in  rigid  body  dynamics.  The  Euler  angles  are 
used  only  to  specify  an  initial  orientation  for  the  body.  The  next  section  contains  formulae 
necessary  to  determine  quarternions  knowing  the  Euler  angles  and  vice  versa. 


5.2.1  Quarternions  in  Terms  of  Euler  Angles 

Let  Gij  be  the  coefficients  of  the  rotation  matrix  G.  The  quarternions  are  then: 
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§4  —  0*5  *  [1  +  (?ii  +  G22  "I"  ^33]  * 

Ql  =  (<^23  ~  Gs2)/^4 
Q2  —  (G31  —  Giz)/Aq^  > 

93  =  (<?12  —  <J2l)/4g4  , 

For  the  special  case  where  =  0,  the  following  logic  is  used: 


(5.21) 

(5.22) 


fi  =  l  +  Gii  :  2  =  1,2,3 

0.5  *  Giz/qi 
0.5  *  G23/92 
0.5  *  (J32/ 93 

end  if 

5.2.2  Euler  Angles  in  Terms  of  Quarternions 


if  (A  /  0)  then 

5i  =  (0.5*A)^  ■  92  =  0.5*^12/91  :  93  = 
else  if  (A  0)  then 

92  =  (0.5  *72)5  :  91  =  0.5*^21/92  :  93  = 
else  if  (A  i=-  0)  then 

93  =  (0.5  *73)^  :  9i  =  0.5*G3i/93  :  92  = 


sin  d  =  2(9294  -  9193) 

There  is  no  loss  of  generality  in  restricting  Q  to  the  range  — 7r/2  <  0  <  7r/2  . 


(5.23) 


tarnl^  = 


tancj)  = 


2(gi92  +  9394) 
2(9l  +  9l)-l 
2(9293  +  9194) 
2(9?  +  94)  -  1 


(5.24) 

(5.25) 


'll)  and  4)  should  be  determined  by  the  ATAN2  function,  to  obtain  values  in  the  correct 
quatrant. 


5.3  The  Full  Equations  of  Rigid  Body  Dynamics 

Equation  (5.8),  (5.9),  and  (5.20)  constitute  a  set  of  10  simultaneous  odes  (ordinary  differential 
equations)  in  time.  The  addition  of  three  more  equations  governing  the  position  vector  of  the 
mass  center, viz. 


fo  =  G’^Uq  (5.26) 

completes  the  set  of  equations  that  are  needed  to  determine  the  trajectory  of  the  body.  These 
equations  are  solved  as  an  initial  value  problem  using  Runge-Kutta  timestepping. 
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5.4  Store  Simulation 


5.4.1  Problem  Statement 

•  Let  [X,Y,Z]  be  an  inertial  frame  (see  Figure  5.1).  This  could,  for  example,  be  fixed  to 
an  aircraft  flying  along  a  steady,  straight  path. 

•  Let  r  be  a  domain  representing  the  farstream,  fixed  to  [X,  F,  Z]. 

•  Let  Ofc  be  a  body  whose  inertial  orientation  and  velocity  with  respect  to  \X,Y,Z]  is 
known.  \X,Y,^  is  a  local  axes  system  attached  to  the  body. 

5.4.2  Solution  Strategy 

First,  it  is  necessary  to  specify  a  set  of  independent  parameters  which  compose  the  solution 
set.  These  are: 

U  The  vector  of  conservation  flow  variables  at  each  grid  point. 

Ugfitnt,  The  grid  velocity  components  at  each  point  in  the  interior  of  the  mesh. 

Xgfiow  The  grid  coordinates  of  interior  nodes. 

q  Quaternions  governing  orientation  of  [X,  F,  Z\. 

uq  =  (uo,  Velocity  components  of  the  mass  center  “0”  of  the  body  in  the 

[X,  F,  Z]  system. 

[Xo,  Fo,  Zq]  —  vq  Coordinates  of  “0”  in  [X,  F,  Z]. 

9?  Rotational  velocity  of  the  body  in  the  [X,  F,  Z]  system. 

fb  Coordinates  of  body  nodes  in  the  [X,  F,  Z]  system. 

We  consider  the  following  two  problems: 

(a)  Given  the  mass  properties  of  Q.b,  compute  its  trajectory  under  the  influence  of  external 
forces  (such  as  gravitational  and  aerodynamic). 

(b)  Given  the  trajectory  f2j„  compute  the  aerodynamic  forces  acting  on  it. 

Next  we  define  the  following  derived  quantities  that  are  important  to  the  solution  methodology. 

^gbody  -  The  grid  velocity  components  of  nodes  on  the  body  in  [X,  F,  Z]  system. 

=  Uq  +  U}b®^b 

Ugbody  •  Grid  velocity  of  nodes  on  the  body  in  [X,  F,  X].  , 

—  G  {q)y'gbody 

Xbody  '•  Function  of  (f6,g,ro) 

=  ro  +  G~’-(g)n 

where  G  is  the  rotation  matrix  relatin  X  to  X  by  X  =  G{q)  X. 

At  each  node  point  in  the  mesh,  the  system  of  equations  that  are  solved  for  the  flow  solver, 
mesh  mover,  and  six-dof  solver  are: 
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Flow  solver: 

ir”  =  F(I7",U;^,C7;^,„,X”)  (5.27) 

Mesh  mover: 

^  gflow 

x” 

gflow 

Six-DOF  solver: 

T 

[uo,Wo,ro,go  =  S{uQ,u}o,ro,qQ,  f^,mext)  (5-29) 

where  vriext  are  the  external  forces  and  moments  respectively.  Equations  (5.27),  (5.28), 
and  (5.29)  constitute  a  set  of  equations,  explicit  in  time,  that  can  be  updated  by  forward  Euler 
or  Runge-Kutta  time  stepping. 

If  we  are  interested  in  solving  only  the  problem  of  specified  trajectory  T(i),  only  equations 
(5.1)  and  (5.2)  need  to  be  solved.  In  addition  we  have: 

55,  *;]’■  =  T(()  (5.30) 

where  t  is  the  current  time  value.  The  body  grid  velocity  and  the  new  position  can  them  be 
computed  from  the  following: 

G-\q^)u^,iody 

r-,+G-\q^)h 

5.5  Node  Movement 

As  the  prescribed  bodies  move  (either  under  captive  or  calculated  trajectories),  the  mesh 
around  the  bodies  must  be  updated  to  avoid  grid  tangling.  This  is  accomplished  in  two  steps: 
1)  grid  node  movement  and  2)  mesh  restructuring.  To  achieve  movement  of  interior  nodal 
points  in  response  to  the  motion  of  the  grid  boundary,  we  appeal  to  an  incompressible  flow 
analogy.  The  advantage  of  this  method  is  that  it  allows  the  nodes  to  flow  through  the  domain 
and  around  the  body  as  it  moved  through  the  field  [16]. 

Assume  that  each  grid  node  represents  a  particle  of  mass  (the  mass  can  be  different  for 
each  node).  We  define  the  ‘density’  of  the  fluid  as  the  mass  of  the  node  divided  by  the  volume 
of  all  connecting  elements.  Initially,  the  density  is  assumed  to  be  unity  (since  each  node’s 


=  (5-28) 
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mass  is  exactly  equal  to  the  volume  of  the  surrounding  elements).  We  then  require  that  the 
fluid  is  incompressible  and  thus  the  density  remains  constant.  Interior  nodes  move  in  response 
to  moving  boundaries  according  to  grid  velocities  calculated  by: 


Vi  =  j3'7pi  (5.31) 

where  /5  is  a  parameter  calculated  to  maintain  overall  stability  of  the  grid  movement  algorithm. 
Given  the  grid  velocities,  actual  node  movement  proceeds  by  integration  of  the  following 
equation: 


(5.32) 


The  node  movement  is  done  using  an  analogy  with  incompressible  flow  that  ensures  that  the 
mesh  moves  consistently  with  the  bodies,  and  that  an  acceptable  mesh  spacing  is  maintained. 
The  node  movement  strategy,  since  it  is  almost  identical  to  the  flow  solver  methodology,  is 
fully  vectorized  in  the  same  manner  as  the  flow  solver. 


5.6  Mesh  Restructuring 

After  several  time  steps,  it  may  be  necessary  to  restructure  the  mesh  by  keeping  the  node 
positions  fixed,  but  changing  the  mesh  connectivity.  We  recover  a  Delaunay  mesh  after  node 
movement  by  performing  the  3D  edge  swapping  procedure  outlined  above.  First  we  update  the 
node  locations  using  the  grid  node-flow  analogy.  Next,  we  visit  groups  of  five  nodes  that  form 
either  two  or  three  tetrahedra  and  swap  the  edges  where  the  Delaunay  criterion  is  violated. 
This  method  has  been  implemented  very  efiiciently  by  keeping  a  hash  table  of  edges  in  the 
mesh  that  have  exactly  three  surrounding  tets  (and  are  candidates  for  swapping  into  two  new 
tets). 
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6  Visualization 


Visualization  is  a  key  piece  of  any  three-dimensional  CFD  simulation,  and  STORESIM  pro¬ 
vides  a  full  suite  of  capabilities  to  visualize  the  time-dependent,  3D  solutions.  The  visualization 
module  of  STORESIM  is  called  STOREVIS.  Users  can  select  any  time  value  to  be  visualized 
(as  long  as  it  was  saved  into  the  restart  file).  If  the  time  value  does  not  exactly  correspond  to 
a  value  in  the  file,  the  solution  values  are  linearly  interpolated  from  the  neighboring  values, 
transparently  to  the  user.  Once  a  time  value  is  chosen,  the  user  can  look  at  (transparent)  pla¬ 
nar  slices,  (transparent)  iso-surface,  particle  traces,  and  boundary  contours  of  various  quantity. 
The  mesh  at  the  given  time  level  can  also  be  viewed. 

The  visualization  package  enables  the  user  to  inspect  the  grid  generation  process  interac¬ 
tively.  It  also  enables  the  user  to  display  selected  regions  for  critical  inspection  and  modifica¬ 
tion  of  the  input  data.  This  is  especially  useful  during  mesh  generation  as  it  provides  the  user 
an  option  to  change  the  data  without  editing  the  input  files. 

The  post-processing  module  enables  the  user  to  display  multiple  slices  of  the  flow  domain. 
It  also  allows  the  user  to  plot  the  selected  variables  on  an  isosurface  of  a  different  variable. 
For  example,  the  user  can  plot  the  variation  of  density  on  an  iso-surface  of  the  pressure.  The 
velocity  vectors  and  particle  traces  option  allow  the  user  to  plot  the  flow  field  and  understand 
the  flow  behavior  in  various  regions  of  the  domain. 
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Results  and  Discussion 


7.1  Preliminary  Remarks 

In  the  previous  sections  we  have  described  important  features  of  the  mesh  generation,  flow 
solver,  six-dof  solver,  and  visualization  modules  of  STORESIM.  In  this  section  we  present  a 
number  of  validation  and  demonstration  problems  which  demonstrate  the  utility  and  flexibility 
of  the  package.  In  each  case  all  of  the  modules  are  used.  In  other  words,  we  start  with 
the  geometry  definition,  then  generate  an  unstructured  mesh  using  TETMESH.  Next,  we 
use  the  integrated  flow  solver  STOREFLOW  (including  the  six-dof  solver  and  mesh  mover 
modules  when  appropriate)  to  obtain  the  flow  solution.  Finally,  the  results  are  post-processed 
using  STOREVIS.  The  problems  have  a  varying  degree  of  complexity  in  terms  of  both  mesh 
generation  and  flow  analysis.  These  problems  have  been  selected  to  test  all  of  the  major 
features  of  the  software. 

7.2  Validation  Examples 

7.2.1  Carter  Problem 

First  we  solve  flow  over  a  flat  plate.  The  free  stream  Mach  number  of  the  flow  is  3.0  and  the 
Reynolds  numbers  is  1000,  and  the  temperature  of  the  plate  is  held  at  IdOOR”.  This  problem  is 
solved  to  validate  the  viscous  flow  solver.  The  problem  is  modeled  using  a  1800  quadrilateral 
elements  with  element  biasing  towards  the  wall.  Figure  7.1  shows  the  [ressure  coefficient  Cp, 
Mach  number,  and  density  variations  along  the  solid  wall  and  in  the  flow  domain.  The  Mach 
contours  show  a  boundary  layer  formation  at  the  isothermal  solid  wall.  Also  this  indicates 
the  flow  causes  a  curved  shock  at  the  leading  edge  of  the  plage.  These  results  agree  with  the 
values  presented  in  [20]. 

Next  we  increase  the  Reynolds  number  to  IE  +  06  and  solved  the  problem  using  the  zero- 
equation  turbulence  model.  Computed  values  of  Cp,  Mach  number  and  density  are  plotted  in 
Figure  7.2.  Comparing  these  with  the  laminar  flow  solution  indicates  that  the  wall  boundary 
layer  for  turbulent  flow,  as  expected,  is  much  thinner  than  the  laminar  boundary  layer. 


7.2.2  Supersonic  Flow  Over  a  Sphere 

As  a  second  test  case,  we  consider  the  problem  of  supersonic  inviscid  flow  over  a  sphere.  The 
radius  of  the  sphere  is  1  unit.  The  radius  of  the  outer  (farfield)  boundary  is  assumed  to  be  10 
units.  The  free  stream  Mach  number  is  assumed  to  be  3.0.  The  inner  and  outer  boundaries  are 
each  represented  by  eight  network  patches  as  shown  in  Figure  7.3.  Each  of  these  16  patches 
has  one  degenerate  edge.  To  produce  a  similar  triangulation  on  each  surface,  we  discretized 
each  edge  using  15  nodes.  This  edge  discretization  produced  188  triangles  on  each  surface 
and  3008  triangles  on  all  16  surfaces.  Surface  triangulation  on  the  inner  boundary  is  shown  in 
Figure  7.4.  Next  we  extruded  the  inner  surface  by  0.25  units  to  produce  5  layers  of  prismatic 
elements.  The  prismatic  elements  are  shown  in  Figure  7.5  The  final  mesh  consists  of  17,156 
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nodes,  80,406  tetrahedral  elements  and  7520  prisms.  The  complete  mesh  generation  process 
required  600  seconds  on  SGI  Indigo  r4000  running  IRIX  5.3. 

Here  we  present  the  results  for  inviscid  flow  analysis.  Mach  number  contours  (see  Figure 
7.6a)  show  the  formation  of  a  bow  shock  ahead  of  the  sphere.  The  standoff  distance  and  the 
stagnation  properties  are  in  good  agreement  with  the  experimental  values.  Velocity  vectors 
(Figure  7.6b)  at  the  midplane  of  the  sphere  in  the  flow  direction  show  a  recirculation  region 
behind  the  sphere.  Pressure  contours  and  density  contours  are  shown  in  Figure  7.6c  and  7.6d 
respectively. 

7.2.3  Double  Ellipsoid 

Next  we  consider  hypersonic  flow  over  double  ellipsoid.  Surface  geometry  of  the  double  ellip¬ 
soid  is  given  in  Figure  7.7.  The  free  stream  Mach  number  for  this  problem  is  8.15  and  the 
angle  of  attack  is  30  deg.  Figure  7.8  shows  the  trimmed  surface  and  their  normals,  edges  of 
the  trimmed  surface,  and  the  edge  discretization.  The  node  distribution  along  each  edge  is 
shown  in  green  circles.  Figure  7.9  shows  the  wireframe  representation  of  the  double  ellipsoid 
surface.  This  flgure  shows  the  cusped  narrow  regions  where  the  top  and  bottom  ellipsoids 
meet.  Such  surfaces  are  especially  challenging  to  automatic  mesh  generators  because  of  the 
surface  curvatures  and  non-unique  normals.  Figures  7.10  and  7.11  show  the  surface  trian¬ 
gulation.  A  total  of  15,794  surface  triangles  were  created  on  both  double  ellipsoid  surface 
and  the  farstream  bounding  box.  After  volume  meshing  the  final  mesh  had  65,608  nodes  and 
392,921  tetrahedral  elements.  The  complete  mesh  generation  process  (including  triangulation 
and  tetrahedralization)  required  5400  seconds  on  SGI  Indigo  r4000  running  IRIX  5.3. 

The  free  stream  Mach  number  and  the  angle  of  attack  are  taken  to  be  equal  to  8.15  and 
30  deg,  respectively.  Inviscid  flow  analysis  results  are  compared  with  those  presented  in  [7]. 
Figure  7.12  shows  the  contours  of  density,  and  coefficient  of  pressure  on  the  double  ellipsoid 
surface.  Also,  plotted  in  Figure  7.12  are  variation  of  density  and  pressure  coefficient  along  the 
symmetry  plane.  Computed  values  are  compared  with  those  presented  in  [7]  and  they  show 
excellent  agreement.  Figure  7.13  shows  the  Mach  contours  at  the  symmetry  plane.  These 
contours  show  a  bow  shock  formation  ahead  of  the  nose  of  the  ellipsoid  and  another  shock 
formation  at  the  intersection  of  the  two  ellipsoids. 

7.2.4  Flow  over  a  3-D  Obstacle 

Here  we  consider  the  problem  of  flow  over  a  three-dimensional  obstacle.  A  Schematic  diagram 
of  the  flow  domain  is  shown  in  Figure  7.14.  Surface  triangulation  is  shown  in  Figure  7.15.  The 
swept  angle  of  the  obstacle  is  30  deg  and  the  plates  are  inclined  at  an  angle  or  15  deg.  The 
free  stream  Mach  number  and  Reynolds  number  are  assumed  to  be  equal  to  10  and  9.0E4-06, 
respectively.  This  flow  is  modeled  as  laminar  nonreacting  flow.  The  mesh  consists  of  10 
138  surface  triangles,  25  619  nodes,  and  146  362  tetrahedral  elements.  The  complete  mesh 
generation  process  required  2200  seconds  on  SGI  Indigo  r4000  running  IRIX  5.3.  Figure  7.16 
shows  the  velocity  at  the  mid  plane.  Figure  7.17  shows  the  Mach  contours. 
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Figure  7.4:  Surface  triangulation  Figure  7,5:  Prismatic  elements  on  one  patch. 
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Figure  7.7:  Schematic  Diagram  of  Double  Ellipsoid  Geometry. 
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Figure  7.8:  Trimmed  surfaces,  surface  normals,  edge  normals,  and  edge  discretization  for  double  ellipsoid  geometry. 


49 


Figure  7.9:  Wireframe  representation  of  the  double  ellipsoid  surface.  This  is  used  as  the  input  and  the  surface 
triangles  are  created  such  that  the  nodes  fall  on  this  surface. 
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Figure  7.10:  Triangulation  of  the  surface  of  the  double  ellipsoid.  Total  number  (both  on  outer  boundary  and  the 
double  ellipsoid  surface)  of  triangles  created  is  equal  to  15,794. 
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Figure  7.11:  Another  view  of  surface  triangles.  This  shows  the  triangulation  on  part  of  the  far-field  boundary. 
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Figure  7.13:  Hypersonic  flow  over  double  ellipsoid.  Mach  contours  at  the  midplane 
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Figure  7.14:  Schematic  diagram  of  three-  dimensionalobstacle. 
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compression  comer  flow  problem. 


7.2.5  Flow  of  Real  Gas  over  Cylinder 

To  demonstrate  the  real  gas  modeling  capabilities  of  STOREFLOW,  we  consider  hypersonic 
flow  over  a  cylinder.  The  real  gas  simulation  modeling  in  STOREFLOW,  instead  of  using 
tables  or  correlations  as  used  by  many  researchers,  employs  exact  thermodynamical  relations. 
The  real  gas  model  assumes  equilibrium  chemistry  and  a  flve  species  model  for  air  modeled  by 
statistical  thermodynamics  [1].  These  five  species  are;  molecular  nitrogen  molecular  oxygen, 
nitric  oxide,  atomic  oxygen  and  atomic  nitrogen. 

The  free  stream  Mach  number  of  the  fluid  is  10.  This  is  a  two-dimensional  test  case  and 
the  finite  element  mesh  shown  in  Figure  7.18  consists  of  2940  quadratic  elements  with  proper 
biasing  close  to  the  cylinder.  In  Figure  7.19  we  compare  temperature  and  Mach  contours 
for  real  gas  flow  analysis  with  those  of  perfect  gas  flow  at  the  same  Mach  number.  The 
temperatures  at  the  cylinder  wall  and  elsewhere  in  the  domain  are  higher  in  the  case  of 
perfect  gas  simulation.  Also,  the  real  gas  case  shows  a  thinner  shock  layer  and  is  closer  to  the 
cylinder.  Contours  the  five  species  are  shown  plotted  in  Figure  7.20. 

7.2.6  Generic  Submunitions  on  CBU 

In  this  case  we  generate  tetrahedral  mesh  for  a  Cluster  Bomb  Unit  (CBU)  with  four  submuni¬ 
tions.  The  finite  element  mesh  consists  of  18  128  surface  triangles  on  the  boundaries.  The  fins 
of  the  CBU  are  modeled  as  zero-thickness  surfaces.  For  such  surfaces,  the  mesh  generator  au¬ 
tomatically  creates  matching  surface,  domain  boundary  and  corresponding  surface  triangles. 
Figure  7.21  shows  surface  triangulation  of  the  CBU.  After  tetrahedralization,  the  mesh  has  25 
795  nodes  and  136  474  tets.  The  complete  mesh  generation  process  required  4200  seconds  on 
SGI  Indigo  r4000  running  IRIX  5.3.  The  difference  in  the  CPU  time  from  previous  problem 
is  due  to  the  complexity  of  the  surface  geometry  in  the  problem. 

Here  we  solved  the  problems  using  the  Euler  module  of  the  flow  solver.  The  freestream 
Mach  number  is  0.85.  Contours  of  density,  pressure,  and  Mach  number  on  the  surface  of  CBU 
are  shown  in  Figure  7.22.  Velocity  vectors  on  the  midplane  are  shown  in  Figure  7.23.  These 
show  recirculation  in  the  cavities. 

7.2.7  F-16 

The  F16  geometry  was  prepared  by  WL  personnel  and  delivered  to  COMCO  (Figure  7.24). 
This  geometry  consists  of  a  detailed  representation  of  the  F16  with  a  fuel  tank  and  a  store 
mounted  under  each  wing.  A  full  set  of  experimental  data  has  been  obtained  for  the  store 
for  various  free  stream  conditions  [12].  The  geometry  was  discretized  using  TETMESH.  The 
mesh  consists  of  28  661  surface  triangles,  119  829  nodes,  719  881  tetrahedral  elements,  and 
95  300  prismatic  elements  surrounding  the  F16,  fuel  tank,  and  the  store.  CPU  time  required 
for  mesh  generation  was  8500  seconds  on  an  SGI  Indigo  r8000  running  IRIX  6.0.  Figure  7.25 
shows  triangular  mesh  on  the  FI  6  surface. 

Figures  7.26  show  contours  of  pressure  coefficient  Cp  on  F-16  and  store  for  an  inviscid  run. 
Here  the  free-stream  Mach  number  is  0.95  and  the  angle  of  attack  is  zero  degrees.  The  solution 
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Figure  7.18  Finite  element  mesh  for  flow  over  a  cyliner  (M  =  10) 
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imire  7.20:  Real  gas  flow  over  a  cylinder  (M  =  10).  Contours  of  Mach  mimber  and  five  species  of  air 
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Figure  7.21:  Surface  triangulation  of  CBU  geometry  with  4  generic  submunitions. 
Fins  of  CBU  are  defined  as  zero-thickness  surfaces. 
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Figure  7.24:  Surface  geometry  of  F-16  with  fuel  tanks  and  generic  store. 
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Figure  7.25:  Surface  triangulation  on  F-16  with  fuel  tanks  and  generic  store. 


was  said  to  be  converged  when  the  non-dimensional  solution  norm  was  reduce  by  4  orders  of 
magnitude.  It  took  6500  time  steps  for  to  obtain  a  converged  solution.  Figure  7.27  shows 
the  comparison  pressure  coefficient  on  the  aircraft  body  at  two  different  free-stream  Mach 
numbers  0.95  and  1.1.  In  Figure  7.28  we  compare  the  Cp  on  the  fins  of  the  store  at  M  =  0.95 
and  1.1.  The  computed  results  predict  the  overall  trend  and  capture  the  flow  behavior.  The 
differences  are  attributed  to  the  fact  that  the  numberical  solutions  were  obtained  using  inviscid 
flow  assumption. 

7.3  Demonstration  Cases 

7.3.1  F-16  with  free  falling  generic  store 

Starting  with  the  converged  solution  from  previous  inviscid  F-16  run  as  the  intial  flow  field, 
we  released  the  store.  The  store  separation  calculations  were  performed  using  the  integrated 
flow-,  six-dof-,  and  mesh-moving,  algorithm  implemented  in  STORESIM.  The  analysis  was 
run  for  another  3250  time  steps.  For  each  time  step  the  location  of  the  store  was  updated 
along  with  any  required  restructuring.  The  initial  conditions  of  the  store  at  the  start  of  the 
separation  are; 

Step:  0  iDU.DVi/dt  =  9.42779 

[get_sixdof .motion] : :  motion  information  at  time:  0 
com.velocity .tilde :  0  0  -0.203 
omega.tilde:  000 
quaternions :  0  0  0  1 
center.of .mass :  339.328  120  63.3017 

The  location  of  the  center  of  mass  and  its  velocity  at  different  time  steps  is  given  in  below. 
The  data  files  and  the  solution  save  files  are  provided  in  the  accompanying  tape. 

Step:  100  iDU.DVi/dt  =  4.28933 

Cget.sixdof . motion] : :  motion  information  at  time:  0.00642502 
com.velocity.tilde :  0.00283813  -0.00048171  -0.266419 

omega.tilde:  7.84291e-05  -7.99302e-05  -5.34408e-05 

quaternions:  1.25381e-07  -1.27756e-07  -8.53401e-08  1 

center.of .mass :  339.328  120  63.3002 

xode(l):  0.00642502 

xode(2:4):  0.00283813  -0.00048171  -0.266419 
xode(5:7):  7.84291e-05  -7.99302e-05  -5.34408e-05 

xode(8:ll):  1.25394e-07  -1.27769e-07  -8.53487e-08  1 

xode(12:14):  339.328  120  63.3002 

Step:  500  iDU.DVi/dt  =  0.28565 

[get.sixdof. motion] : :  motion  information  at  time:  0.0319988 
com.velocity.tilde:  0.0141818  -0.00259118  -0.518533 
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Figure  ..26:  F16p(3Dontour. 


Figure  7.27:  Comparison  of  computed  values  of  Cp  over  F-16  with  experimental  data 
at  M  =  0.95  and  1.1 
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omega_tilde:  0.000394509  -0.000402629  -0.000271088 

quaternions:  3.14511e-06  -3.20707e-06  -2.15391e-06  1 

center_of_mass:  339.328  120  63.2901 

xode(l):  0.0319987 

xode(2:4):  0.0141818  -0.00259118  -0.518533 
xode(5:7):  0.000394509  -0.000402629  -0.000271088 

xode(8:ll):  3.14512e-06  -3.20709e-06  -2.15392e-06  1 

xode(12:14):  339.328  120  63.2901 

Step:  1000  iDU.DVi/dt  =  0.194849 

[get _sixdof .motion] : :  motion  information  at  time:  0.063968 
com. velocity .tilde:  0.0283269  -0.00559457  -0.833152 

omega_tilde:  0.000792428  -0.000813205  -0.000550542 
quaternions:  1.26306e-05  -1.2915e-05  -8.71066e-06  1 
center.of.mass:  339.328  120  63.2685 
xode(l):  0.063968 

xode(2:4):  0.0283269  -0.00559457  -0.833152 
xode(5:7):  0.000792428  -0.000813205  -0.000550542 
xode(8:ll):  1.26306e-05  -1.2915e-05  -8.71067e-06  1 
xode(12:14):  339.328  120  63.2685 

Step:  1500  iDU.DVi/dt  =  0.154484 

[get.sixdof. motion] : :  motion  information  at  time:  0.0959398 
com. velocity .tilde:  0.0424375  -0.00895945  -1.14731 

omega.tilde:  0.001189  -0.00122945  -0.00083683 
quaternions:  2.84727e-05  -2.92358e-05  -1.9792e-05  1 
center.of.mass:  339.328  120  63.2369 
xode(l):  0.0959398 

xode(2:4):  0.0424375  -0.00895945  -1.14731 
xode(5:7):  0.001189  -0.00122945  -0.00083683 
xode(8:ll):  2.84727e-05  -2.92358e-05  -1.9792e-05  1 
xode(12:14):  339.328  120  63.2369 

Step:  2000  iDU.DVi/dt  =  0.127969 

[get.sixdof. motion] : :  motion  information  at  time:  0.127914 
com.velocity.tilde :  0.0565076  -0.0126205  -1.4611 

omega.tilde:  0.00158002  -0.00164992  -0.00112796 
quaternions:  5.06172e-05  -5.2247e-05  -3.54926e-05  1 
center.of.mass:  339.328  120  63.1952 
xode(l):  0.127914 

xode(2:4):  0.0565076  -0.0126205  -1.4611 
xode(5:7):  0.00158002  -0.00164992  -0.00112796 
xode(8:ll):  5.06172e-05  -5.2247e-05  -3.54926e-05  1 
xode(12:14):  339.328  120  63.1952 
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step:  2500  iDU.DVi/dt  =  0.106743 

[get_sixdof .motion] : :  motion  information  at  time:  0.15989 


com_velocity_tilde : 
omega.tilde : 
quaternions : 
center.of .mass : 
xode(l) : 
xode(2:4) : 
xode(5:7) : 
xode (8:11) : 
xode(12;14) : 


0.0705241  -0.0165259  -1.77459 
0.00196262  -0.00207458  -0.00142251 
7.89492e-05  -8.20148e-05  -5.58772e-05 

339.328  120  63.1434 
0.159892 

0.0705241  -0.0165259  -1.77459 
0.00196262  -0.00207458  -0.00142251 
7.89492e-05  -8.20148e-05  -5.58772e-05 

339.328  120  63.1434 


1 


1 


Step:  3000  iDU.DVi/dt  =  0.0918054 

[get.s ixdof. motion] : :  motion  information  at  time:  0.191867 


com.velo city. tilde : 
omega.tilde : 
quaternions : 
center.of .mass : 
xode(l) : 
xode (2; 4) : 
xode (5: 7) : 
xode (8: 11) : 
xode (12: 14) : 


0.0844736  -0.020651  -2.08782 
0.00233644  -0.0025036  -0.00172001 
0.000113327  -0.000118608  -8.0995e-05  1 

339.328  120  63.0817 
0.19187 

0.0844736  -0.020651  -2.08782 
0.00233644  -0.0025036  -0.00172001 
0.000113327  -0.000118608  -8.0995e-05  1 

339.328  120  63.0817 


Step:  3250  iDU.DVi/dt  =  0.086224 

[get.s ixdof. motion] : :  motion  information  at  time:  0.207855 
com. velocity. tilde:  0.0914193  -0.0227924  -2.24434 

omega.tilde:  0.00252038  -0.00271971  -0.0018699 
quaternions:  0.000132742  -0.000139485  -9.53437e-05  1 
center.of .mass:  339.328  120  63.0471 
xode(l):  0.207859 

xode (2: 4):  0.0914193  -0.0227924  -2.24434 
xode (5: 7):  0.00252038  -0.00271971  -0.0018699 
xode (8: 11):  0.000132742  -0.000139485  -9.53437e-05  1 
xode (12: 14):  339.328  120  63.0471 


The  displacement  of  the  store  at  the  end  of  3250  time  steps,  after  converting  the  units  to 
the  scaled  down  model  units,  is  the  vertical  displacement  of  the  store  is  11  inches. 

in  [26] 
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7.3.2  Hypersonic  2-Stage  Case 

Details  of  this  demonstration  case  are  presented  in  [6] 
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8  Conclusions 


The  primary  goal  of  this  effort  was  to  couple  these  disciplines  into  a  unified  computational 
environment  for  analyzing  dynamic  separation  of  weapons,  stages,  and  crew  escape  vehicles 
from  supersonic  and  hypersonic  parent  vehicles.  The  major  outcome  of  this  research  and 
development  effort  is  a  software  package  called  STORESIM.  The  STORESIM  package  is  a 
tightly  integrated  group  of  modules  (TETMESH,  STOREFLOW,  STOREVTS)  that  performs 
grid  generation,  static  fiow  simulation,  six  degree-of-freedom  rigid  body  dynamics,  dynamic 
flow  analysis,  body  movement,  grid-node  movement,  grid  restructuring,  and  time-dependent 
results  visualization. 

Key  accomplishments  of  this  effort  are: 

1.  Mesh  Generation:  A  grid  generator  is  essentially  useless  without  a  powerful  technique  for 
defining  the  geometry  to  be  discretized.  To  this  end  we  have  developed  an  unstructured 
mesh  generation  module,  TETMESH,  that  is  capable  of  reading  the  geometry  in  a  very 
simple  network/trimmed  surface  format.  Salient  features  of  this  module  include: 

•  Triangulation:  Surface  triangulation  is  based  on  the  Dulanay  criterion.  The  surface 
triangulator  is  capable  of  handling  complex  curved  surfaces  with  cusped  edges. 

•  Volume  Gridding:  The  volume  mesher  generates  tetrahedral  elements  with  an  op¬ 
tion  of  generating  prismatic  elements  extruded  from  the  surfaces.  The  tetrahedral- 
ization  is  also  based  on  Dulanay  criterion  with  proper  face  enforcement  to  fill  the 
voids. 

•  Zero-thickness  Surface:  Thin  fins  protruding  from  a  missile  can  be  approximated 
as  single  surfaces  during  data  preparation.  TETMESH  automatically  generates  a 
matching  surface  to  create  a  barrier  and  prevent  fluid  movement  across  this  surface. 

•  Adaptivity:  The  meshing  module  is  capable  of  refining/unrefining  the  mesh.  Re¬ 
gions  of  maximum  solution  gradients  are  automatically  refined  to  better  capture 
solution  in  critical  regions. 

2.  Flow  Solver:  The  STORESIM  package  has  a  compressible  fluid  flow  solver  module  called 
STOREFLOW.  Key  features  of  the  STOREFLOW  module  are: 

•  An  option  for  ALE  (Arbitrary  Lagrangian  Eulerian )  computations.  This  essentially 
makes  it  possible  to  use  moving  meshes  to  model  multiple  moving  bodies.  This 
new  ALE  formulation  is  extremely  robust,  and  can  even  produce  accurate  results 
for  bodies  launched  at  high  speed  into  ambient  air. 

•  A  very  sophisticated  far  stream  boundary  condition  designed  to  prevent  spurious 
energy  build  up  and  consequent  reflections  at  transonic  speeds.  This  approach  also 
includes  modifications  to  account  for  moving  mesh  boundaries. 

•  Euler  and  Navier-Stokes  versions.  Turbulence  modeling  based  on  a  simple  zero 
equation  model  has  been  tested  for  two-dimensional  flows.  The  extension  to  three- 
dimensional  turbulence  modeling,  although  not  implements  should  be  straightfor¬ 
ward. 
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•  Real  gas  capability,  based  on  a  5  species  (Nitrogen,  Oxygen,  Atomic  Nitrogen, 
Atomic  Oxygen  and  Nitric  Oxide)  model  for  dissociated  air. 

•  Multiple  element  types  (tetrahedra,  Hexahedra  or  Prisms).  Additional  element 
types  can  be  easily  incorporated  into  the  code. 

3.  Six-DOF  Solver:  STOREFLOW  is  fully  integrated  with  a  six  degree-of-freedom  (6  dof) 
solver  for  rigid  body  dynamics,  and  a  meshmover,  that  moves  the  mesh  in  response  to 
body  motions. 

4.  Node  Movement:  An  innovative  grid  movement /restructuring  method  based  on  the  grid- 
node-flow  analogy  and  generalized  edge  swapping  in  three-dimensions.  The  node  move¬ 
ment/restructuring  is  required  when  modeling  body  movements  in  the  flow  field. 

5.  Visualization:  The  visualization  package,  STOREVIS,  is  capable  of  post-processing 
structured/unstructured  grids  and  producing  color  plots  of  slices,  isosurfaces,  velocity 
vectors,  and  particle  traces. 

8.1  Future  Work 

STORESIM  has  been  used  to  solve  several  complex  problems.  In  every  case,  all  relevant 
modules  of  the  STORESIM  package  are  used.  During  the  development  and  final  validation 
process,  each  module  has  been  tested  separately  and  as  a  complete  package  for  accuracy  and 
robustness.  Based  on  these  tests,  we  strongly  feel  that  STORESIM  is  capable  of  generating 
unstructured  meshes  on  complex  geometries  with  very  little  user  intervention.  The  flow  solver 
and  six-dof  solvers  are  very  stable  and  produce  accurate  results.  Since  this  was  a  major 
development  effort  aimed  at  alleviating  the  problems  with  mesh  generation  and  modeling 
multi-body  CFD  problems,  entire  effort  was  put  in  implementation  of  accurate  and  robust 
techniques.  However,  there  are  several  areas,  such  as  interfacing  with  third  party  flow  solvers, 
implementation  of  better  turbulence  models,  that  need  additional  work. 
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A  Appendix  -  A 

A.l  Flux  Jacobian  for  the  Navier  -  Stokes  Equations 

The  Navier-Stokes  equations  expressed  in  terms  of  conservation  variables  are: 

+  =  (A.1) 

Ff  is  the  Euler  Flux  Vector.  (A.l)  can  be  re-written  as: 


u^t  + 


(A.2) 


where  Ai  =  F^tt  =  Euler  Flux  Jacobian.  Formulation  for  Ai  can  be  found  in  many  standard 

i,U 

references. 


A.2  Flux  Jacobians  in  Terms  of  Entropy  Variables 

The  Entropy  variables  are  defined  as: 


=  (ps)u  (A.3) 

where,  p  is  the  density  of  the  fluid,  and  s  is  the  entropy.  Use  of  these  variables  symmetrizes 
the  Navier-Stokes  Equations: 

(A.l)  becomes 


where 


A,Vt  +  AiVi=[kijVi,]. 

L  J  ^2 


Aq  -  U^y 

Ai  =  Ffy  =  FfyU^y  =  AjAo 

Ff  =  kijVij 


(A.4) 


(A.5) 


The  matrices  AQ,Ai,kij  are  defined  in  term  of  the  entropy  (or  conservation  )  variable  in 

[20]. 

A.3  Streamline  Operator  Isupg 

Isupc  =  /  W,i)  .  [r\{At^^  Vi}dn  (A.6) 
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where 


A.4 


where, 


(Ai-uf  I)  =At^^-Ao 

Vector  of  entropy  variables  =  (— p  s)jj 
U,v 

lA-^\Af^^Af^^\-y^ 

Shock  Operator  Ishock 

^shock  ~  f  ^  ■  {^d  ~~  0/U^A(^^idQi 
J  f2 


= 


V-r  = 


ys)  ■  _ 

(Af^E  Vi)  ■  [r]  A-\Af^^  V,i 
(Vk-AoV,,,) 


=  scalar 
=  scalar 


B  Appendix  -  B 

B.l  Source  Code  Distribution  and  Executable  Build  Environment 

Source  code  with  sample  input  files  and  user  manual  has  been  delivered  to  WPAFB  and 

WPAFB  reserves  the  rights  to  distribute  the  software  developed  under  this  contract. 

B.2  Building  STORESIM 

1.  Change  directory  to  scripts  ‘$  cd  scripts’ 

•  Look  at  the  file  sysdep.  This  file  contains  all  system  dependent  information.  Care¬ 
fully  go  over  the  various  compilers,  options  and  paths  specified  here  and  make 
corresponding  modifications  for  your  system.  At  the  minimum,  you  will  need  to 
modify  the  variable  ROOT  which  should  be  set  to  the  pathname  for  storesim  di¬ 
rectory  (i.e.,  the  path  that  you  get  by  doing  ’pwd’  in  storesim  directory). 

•  Check  the  path  specified  for  PERL.  If  it  needs  to  be  modified,  then  make  the  change 
in  sysdep  and  also  make  similar  change  in  file  w4  (the  first  line). 

2.  All  the  makefiles  in  storesim  are  based  on  GNU’s  make.  Therefore,  a  copy  of  this  make  is 
also  included  in  the  scripts  directory  by  the  name  of  makestore_xxx.  You  must  include 
the  path  to  this  executable  in  your  shell  variable  path  (in  your  .cshrc/.tcshrc  file). 

3.  To  make  the  code,  go  to  the  storesim  directory  and  just  say  makestore_xxx. 

B.3  Running  STORESIM 

Please  refer  to  the  accompanying  User  Manual. 
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C  Appendix  -  C 

C.l  User  Manual 
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COMCO 


STORESIM/TETMESH  User  Manual 


1.0  Introduction: 

This  manual  explains  the  basic  steps  necessaiy  to  generate  tetrahedral  meshes  using 
TETMESH  and  perform  the  flow  analysis  using  STORESIM  TETMESH  can  also  be  run  from  the  main 
menu  of  STORESIM.  Please  note  that  this  is  a  preliminary  version  of  the  manual  and  does  not  explain  all 
the  features.  A  complete  description  of  all  features  and  procedures  to  select  various  options  will  be  made 
available  in  future  releases. 

In  the  following,  several  commands  have  been  described  that  invoke  various  modules.  In  the 
definition  of  those  commands  the  variable  $STORESIM_DIR  needs  to  be  set  using  die  standard  setenv 
command  and  its  value  is  the  path  name  to  the  storesim  (root)  directory.  For  example  to  run  STORESIM, 
the  instruction  would  be  to  execute  $STORESIM_DIR/storesim. 


2.0  Input  Geometry  Definition: 

TETMESH  accepts  geometry  information  in  the  form  of  both  IGES  trimmed  surface  format 
and  TETMESH  E3S  trimmed  surface  format.  In  addition,  TETMESH  is  capable  of  converting  Network  (or 
panel  code)  format  to  TETMESH  E3S  trimmed  surface  format. 

2.1  Procedure  to  Convert  a  Network  File  to  ESS  File: 

In  this  section  we  explain  the  steps  involved  in  converting  a  Network  (or  panel  code)  format 
file  to  TETMESH  E3S  trimmed  surfece  format.  First,  TETMESH  processes  data  from  a  Network  file  (say 
model.net)  and  creates  a  file  netlels.  The  data  in  this  file  is  in  E2S  format  and  needs  to  be  converted  to 
E3S  format,  e2s-to-e3s.pl  processes  the  contents  in  net2e2s  file  and  outputs  the  data  in  E3S  format.  The 
necessary  commands  to  convert  the  data  from  Network  format  to  E3S  format  are  given  below. 

a  Tnr.r>nvffrtthftmr)del.nettonet2e2sfonnat  tvue:  $STORESIM_DIR/tet -N  model.net 

h  Tn  convert  net2e2s  to  mnde1.e3s  type:  $STORESIM_DIR/e2s-to-e3s.pl  <  net2e2s  >  modeLeSs 

2.2  Sample  Network  file: 

This  file  contains  the  surface  information  in  the  form  of  a  set  of  (m  x  n)  xyz  coordinates  for 
each  surface.  This  does  not  allow  for  trimmed  surfaces.  Description  of  a  sample  ’Net’  file  is  given  below: 

<Title> 

Number  of  Surfaces 

<Nrunber  of  surfaces> 

<n  >  <n  >  <surface> 
i-l  j-l 

x(l,l)y(l,l)z(l,l) 
x(n._j,l)  y(n._j,l)z(n._j,l) 
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Next  Surface 


Sample  Network  file: 


Sphere  within  a  sphere 
Number  of  Surfaces 
16 

10  10  inner-0 

1.000000  0.000000  0.000000 
0.984808  0.173648  0.000000 

0.766044  0.642788  0.000000 
0.642788  0.766044  0.000000 

10  10  inner-1 

10  10  inner-16 
0.000000  1.000000  0.000000 
-0,173648  0.984808  0.000000 
-0.342020  0.939693  0.000000 


2.3  Sample  £38  file: 

The  E3S  format  is  an  internal  format.  As  explained  in  the  previous  section  ,  the  Network 
files  can  be  converted  into  E3S  format  to  perform  the  mesh  generation  and  flow  analysis  .  In  E3S  format, 
the  geometry  is  defined  by  support-surfaces,  edges,  and  trimmed-surfaces.  The  support-surfaces  are 
defined  by  wireframes  of  (m  x  n)  points.  A  set  of  xyz-points  defining  a  triimning  curve  denotes  an  edge.  A 
trinraied-surface  is  defined  by  a  set  of  boundaries  (a  boundary  is  a  set  of  edges).  The  flow  boundary 
conditions,  component  definition,  material  properties  and  component  trajectory  information  is  appended  to 
complete  the  E3S  format  file.  A  sample  E3S  file  is  given  below. 

/*  support  surfaces  */ 

{SUPPORT-SURFACE  ssl 
{BILINEAR  10  10 
{ 

10  0 

0.984808  0.173648  0 
0  0  1 

} 

} 

} 

/*  boundary  segments  */ 

{TRIMMED-SURFACE  inner-0  ssl  +OUTWARD-NORMAL 
{ 

{  +el  -e2  +e3  ) 

} 


/*  intersection  segments  */ 

{EDGE  el  CURVED  inner-0  inner-1 
{NODES  5  {DISTRIBUTION  INPUT}} 
{ 

0  10 
0  0  1 
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} 

} 

/*  walls  */ 

{BOUNDARY-CONDITION  bcl  SOLID-WALL-EULER 

{0.00027648  1001}  /*  (Temperature,  ul,  u2,  u3,  Density)  */ 

{inner-0  inner-1  inner-2  inner-3  inner-4  inner-5  inner-6  inner-7} 

} 

{COMPONENT  s6dof  SIXDOF-SOLVER-CALCULATE 

{inner-0  inner-1  inner-2  inner-3  inner-4  inner-5  inner-6  inner-7} 

} 

# 

#  The  following  PROGRAM  pi  defines  the  initial  conditions  and  mass/moment-of-inertia 

#  properties  for  the  COMPONENT  s6dof  (this  can  also  be  used  as  a  template) . 

# 

PROGRAM  pi  [ 
s6dof_accg  =  9.81; 
s6dof_amass  =  1.0; 
s6dof_aixx  =  1.0; 
s6dof_aiyy  =  1.0; 
s6dof_aizz  =  1.0; 
s6dof_aixy  =  0.0; 
s6dof_aixz  =  0.0; 
s6dof_aiyz  =  0.0; 

proc  s6doLpiotion  0  { 
time  =S1; 
psi_dot  =  0.0; 
theta_dot  =  0.0; 
phi_dot  =  0.0; 
s6do£_xc  =  0; 

s6dof_yc  =  0; 

s6dof_zc  =  0  ; 

s6dof_vxc  =  5.0; 

s6dof_vyc  =0.0; 

s6dof_vzc  =0.0  ; 

# 

#  YAW  and  YAW-RATE 

# 

s6dof_psi  =  0.0; 

s6dof_psi_dot  =0.0; 

# 

#  PITCH  and  PITCH-RATE 

# 

s6dof_theta  =  0.0; 

s6dof_theta_dot  =  0.0; 

# 

#  ROLL  and  ROLL-RATE 

# 

s6dof_phi  =  0.0; 

s6dof_phi_dot  =  0.0; 

}; 

3.0  Graphics  Server  Setup: 

The  graphics  part  of  STORESIM/TETMESH  is  based  on  the  Motif  and  XI 1  standards.  A 
proper  setup  of  your  workstation  requires: 

1.  Setting  the  environment  variables  DISPLAY  to  the  name  of  your  terminal  by  setenv 
DISPLAY  your_term:0 
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2.  Providing  a  set  of  defeult  bindings  in  the  .Xdefaults  file.  Apart  from  selecting  the  fonts, 
colors  and  locations  of  various  widgets,  the  .Xdefaults  file  also  enables  the  user  to  move  the  objects 
displayed  in  the  viewport  Suggested  defeults  for  the  TETMESH  are  given  in  Table  1.  In  most  cases,  it  is 
suitable  to  copy  the  provided  .Xdefaults  file  to  your  home  directory  under  another  file  name,  and  then 
append  this  file  to  your  existing  .Xdefaults  file  to  preserve  the  contents  of  the  existing  .Xdefaults  file.  These 
additional  .Xdefaults  specification  are  also  repeated  in  a  companion  files  $STORESIM_DIRyXdefaults-tet. 

3.1  Mouse  Button  Functions 

The  user  may  position  the  objects  displayed  in  the  viewport  window  by  pressing  one  of  the 
three  mouse  buttons.  The  definitions  of  these  buttons  corresponding  to  the  .Xdefaults  setting  shown  in 
Table  1  are: 


1.  Button  1  (Left)  Pan:  Move  (drag)  objects  in  the  window. 

2.  Button  2  (Middle)  Trackball:  This  mode  allows  the  objects  to  be  rotated  in  the  window. 

3.  Button  3  (Right)  Zoom:  Moving  the  mouse  horizontally  performs  the  zoom  action, 
basically  enlarging/reducing  the  size  of  the  picture  within  the  window. 


These  defaults  are  for  stand  alone  tet  package 
form  names  are: 

output_form,  exit__form,  edge_e<iitor 


TETMESK*background :  sgi teal 

TETMESH*selectColor :  red 

TETMESH *shadowThickness :  2 

TETMESH*highlightTliickness :  0 

TETMSSH*imagetet. translations ;#override  \n\ 
<BtnlDown>:  Pan(d)  \n\ 

<BtnlUp>:  Pan(u)  \n\ 

<BtnlMotion>:  Pan(m)  \n\ 

<Btn2Down>;  Trackball (d)  \n\ 

<Btn2Up>:  Trackball (u)  \n\ 

<Btn2Motion>:  Trackball (m)  \n\ 

<Btn3Down>:  Zoom(d)  \n\ 

<Btn3Up>:  Zoom{u)  \n\ 

<Btn3Motion>;  Zoom(m)  \n\ 

!  Tetmesh  control  form  and  viewport 
TETMESH *tetview. def aultPosition:  false 

TETMESH  *  tetview . x :  650 

TETMESH  *  t  e t vi ew . y :  550 

TETMESH  *  t  e tvi ew . wi dth :  500 

TETMESH*tetview. height :  500 


TETMESH* tet_form . def aultPosition :  false 

TETME  S  H  *  t  e  t_f  o  rm .  X :  0 

TETMESH*  tet_form.y:  0 
TETMESH*tet_form.  width:  450 

Table  1:  Contents  of  .Xdefaults  file  for. 


4.0  STORESIM  Solver  Options: 

STORESIM  is  a  multibody  CFD  flow  solver  capable  of  handling  moving  boundaries  and 
rigid  body  motions.  The  Move  Components  Form  provides  a  list  of  valid  logical  combinations  of  solvers 
that  can  be  used  in  the  simulation. 
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5.0  Performing  a  Simulation: 

In  this  section  we  give  a  quick  step-by-step  procedure  to  run  STORESIM/TETMESH.  As 
mentioned  in  the  Introduction,  we  do  not  try  to  explain  all  existing  features  and  capabilities  of  the 
STORESIM/TETMESH  program.  However,  for  this  beta  release,  the  steps  described  will  enable  the  user  to 
generate  the  mesh  from  an  E3S  file,  run  the  flow  solver  and  postprocess  the  solution. 

To  start  the  execution  type:  storesim  and  hit  <retum/enter>. 

This  brings  up  die  main  menu  with  the  following  options: 

Read  Files:  To  read  in  restart  or  save  files. 

Viewport :  To  postprocess  the  solution. 

Run  Storeflow:  Execute  the  STORESIM  flow  solver. 

Run  Tetmesh:  Execute  TETMESH  to  generate  mesh  and  execute  STORESIM  flow  solver. 

Save:  To  save  the  solution  to  a  restart  file. 

Done:  End  the  solver  session. 

5.1  Read  Files: 

The  Read  Files  Control  form  consists  of  the  Standard  Motif  Selection  form.  This  is  used  to 
select  and  read  the  solution  file  for  postprocessing.  The  Read  Files  Control  form  is  provided  with  options  to 
read  the  files  in  Tetmesh  grid  file  format,  Storeflow  Solution  File  format,  COBALT  .pix  file  format,  and 
COMCO  restart  file  format.The  form  consists  of  a  scrollable  list  of  file  names,  directories,  input  text  field 
and  Select  and  Filter  buttons. 

5.2  Viewport: 

After  loading  the  solution  file,  by  clicking  on  the  Viewport  button  a  postprocessing  window 
appears  on  the  screen.  This  window  has  several  options  for  postprocessing  the  results.  The  Control  button 
on  the  top  of  this  window  brings  up  the  post-processor’s  primary  control  panel.  This  provides  common 
options,  such  as  selecting  the  background  color,  mesh  display,  color  bar  display,  3D  projection  type,  etc. 
Many  of  these  options  are  activated  by  toggle  switches.  By  clicking  on  the  Load  Results  button  different 
components  can  be  activated  for  postprocessing.  To  activate  (select)  a  component,  click  on  the  name  the  of 
the  component.  The  buttons  act  as  toggle  switches  and  can  be  also  used  to  deactivate  a  component  Only 
active  components  can  be  displayed  in  the  viewport. 

To  activate  die  slicing  plane  feature,  hit  the  Slice  button  on  die  Viewport  menu  or  the  one  in 
the  Control  Panel.  Using  this  option,  scalar  data  may  be  visualized  by  intersecting  the  grid  with  planes 
which  are  colored  according  to  interpolated  values  of  the  desired  solution  component. 

Whereas  slicing  planes  display  scalar  data  according  to  the  geometric  location,  isosurfaces 
are  defined  by  the  data  itself.  The  loci  of  points  with  the  same  data  values  are  an  arbitrarily  shaped  surface 
in  space,  perhaps  widi  disconnected  regions. 

Velocity  vectors  button  enables  the  user  to  plot  the  velocity  vectors  at  the  nodes.  These 
vectors  are  scaled  relative  to  the  maximum  normalized  velocity  in  the  flow  domain  to  provide  an  indication 
of  the  change  in  magnitude  of  the  velocity  vector  at  various  points  in  the  domain.  The  remaining  buttons 
displayed  at  the  top  of  the  vie-wport  window  are  not  yet  completely  functioning  for  the  beta  release. 
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5.3  Run  Storeflow: 

After  loading  the  geometry  file  and  the  mesh,  by  clicking  on  this  button  the  user  executes 
the  STORESIM  solver.  The  link  to  this  button  is  not  yet  properly  established  for  this  release.  However,  for 
this  release,  the  flow  solver  can  be  accessed  through  die  TETMESH.  This  is  explained  in  the  next  section. 

5.4  Run  Tetmesh: 

Choosing  this  option  from  die  STORESIM  main  menu  accesses  the  TETMESH  gird 
generator.  This  brings  up  a  standard  file  selection  form  to  read  the  geometry  information.  The  file  must  be 
in  the  ESS  format  and  have  .eSs  extension.  Once  the  file  has  been  selected,  TETMESH  preprocesses  the 
data  to  verify  the  format.  After  successfully  reading  the  data,  TETMESH  opens  the  following:  a  Viewport 
window  to  display  the  object,  a  Component  Parameters  widget,  and  a  TETMESH  menu  widget.  In  this 
section  we  eiqilain  the  options  and  describe  die  steps  involved  in  generating  the  grid  and  running  the  flow 
solver. 

5.4.1  Viewport: 

The  TETMESH  viewport  is  different  from  the  STORESIM  viewport  described  in  Section 

5.2  The  purpose  of  the  TETEMSH  viewport  is  to  display  the  object  (components)  and  the  mesh.  On  the 
other  hand,  the  STORESIM  viewport  is  designed  to  meet  additional  needs  of  postprocessing  the  solution. 

5.4.2  TETMESH  Menu: 

This  is  die  main  control  form  from  which  the  display  options,  mesh  generation  option  and 
flow  solvers  are  controlled.  At  any  given  time,  all  active  options  are  highlighted  and  the  non-active  options 
are  dimmed  (or  grayed  out).  In  the  control  form  all  options  related  to  a  particular  task  are  grouped  together 
and  are  separated  by  a  dividing  line.  In  the  remainder  of  this  section  we  explain  the  options. 

5.4.2.1  Viewport  Options: 

Surface  Edges:  This  is  a  toggle  used  to  display/clear  the  selected  surfeces. 

Surface  Names:  This  is  a  toggle  used  to  print/erase  the  names  of  selected  surfaces. 

Edge  Names:  This  is  a  toggle  used  to  print/erase  the  names  of  selected  edges. 

Plot  Wireframe:  This  is  a  toggle  used  to  display/erase  the  wireframe  mesh  for  the  surfeces. 

Two-D  Mode:  Display/Erase  the  wireframes  of  selected  surfaces  in  Non-uniform  Parameter  space 
(Parametric  map  of  xyz  in  UV  space). 

Edge  Discretization:  Displays  the  location  of  the  nodes  on  the  selected  surfaces. 

Triangles:  Displays  the  surface  triangles  on  the  selected  surfaces. 

Nodes:  Displays  nodes  on  the  selected  surfaces. 

Non-Unif.  UV-Space:  (Not  Active) 

Unenforced  Faces:  Display  all  enforced  faces. 

Normal:  Display  the  surface  normals  on  selected  surfaces. 

Flip:  Reverse  the  orientations  of  surface  normals  of  selected  surfaces. 

Boundary  Condition:  Print  the  specified  boundary  conditions  on  selected  surfaces. 
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Special  Plot:  For  future  use. 

Tetrahedra  (dynamic^:  Display  all  elements  with  aspect  ratio  less  than  the  minimum  aspect  ratio.  (If  the 
aspect  ratio  of  all  elements  is  greater  than  the  minimum  value,  plot  the  entire  mesh.) 

Tetrahedra  (fasti:  Display  the  tetrahedral  mesh. 

Aspect:  Allowable  minimum  aspect  ratio. 

Reset:  To  reset  the  object  orientation  in  the  viewport. 

S.4.2.2  Surface  Selector: 

Surface  Selector:  This  contains  a  scrollable  list  of  all  selectable  surfaces.  All  surfeces  can  be  selected  by 
clicking  the  left  button  of  the  mouse  on  AU  and  by  cbcking  on  De-sel,  all  surfaces  will  be  de-selected.  All 
selected  components  (surfaces)  will  be  highlighted  and  individual  components  can  be  selected/de-selected 
by  clicking  the  left  button  of  the  mouse  on  the  name  of  the  surfece. 

Surface  Spacing:  Provides  additional  options  to  discretize  the  surface,  (will  be  explained  at  the 
short-course). 


5.4.2.3  Edge  Selector: 

Edges  In  Employ:  This  contains  a  scrollable  list  of  all  selectable  edges.  All  edges  can  be  selected  by 
clicking  the  left  button  of  the  mouse  on  All.  By  clicking  on  De-sel  all  surfaces  will  be  de-selected.  All 
selected  edges  will  be  high-lighted  and  individual  components  can  be  selected/de-selected  by  clicking  the 
left  button  of  the  mouse  on  the  name  of  the  edge. 

Edge  Editor:  Provides  additional  options  to  discretize  the  edges,  (will  be  explained  at  the  short-course). 

5.4.2.4  Triangulate  Surfaces: 

Htt  Everything  fSurfi:  This  is  a  one  step  option  to  triangulate  all  surfaces  in  the  geometry.  This  first 
discretizing  the  edges  and  then  merges  the  nodes  on  matches  edges.  Next  all  surfaces  are  triangulated. 
These  steps  can  be  carried  out  sequentially  by  selecting  the  Discretize  Edges  and  the  options  provided  under 
Triangulate  Surfaces.  The  usages  of  these  options  will  be  explained  at  the  short-course. 

5.4.2.5  Tetrahedralize  RegionsrDo  Everything  fSurfi:  This  is  a  one  step  option  to  tetrahedralize  all 
surfeces  in  the  geometry.  This  first  discretizing  the  edges  and  then  merges  the  nodes  on  matches  edges. 
Next  all  surfaces  are  triangulated.  These  steps  can  be  carried  out  sequentially  by  selecting  the  Discretize 
Edges  and  the  options  provided  under  Triangulate  Surfaces. 

5-4.2.6  SoJyer; 

Move  Components...:  This  enables  the  user  to  select  the  type  of  solver,  time  integration  scheme,  and  other 
solver  parameters.  The  Move  Components  button  in  this  widget  fires  the  appropriate  solver. 

Smooth  Mesh...:  This  widget  provides  the  user  with  different  options  for  selecting  the  mesh  smoothing 
parameters.  In  this  release  only  Node  Flow  algorithm  is  used  for  computing  the  node  velocities  and 
updating  the  mesh. 

Flow  Params...:  The  reference  (farfield/inflow)  flow  parameters  are  displayed  in  this  widget.  Also  some 
run  control  parameters  (CFL  and  number  of  time  steps)  are  q)ecified  here. 

Init  Flowfield:  This  button  is  used  to  initialize  the  flow  field  to  inflow  conditions.  The  existing  solution  (if 
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any)  will  be  lost  and  the  flow  is  initialized  to  inflow  values  on  the  entire  flow  domain. 

Fire  Storeflow:  Only  the  flow  solver  (without  mesh  movement  and  rigid  body  motion)  is  executed  if  the 
user  clicks  on  this  button. 

S.4.2.7  Output  Control; 

This  form  can  be  used  to  save  the  database  into  ESS  format  or  COMCO  restart  format.  The 
COMCO  restart  file  format  can  then  be  used  to  postprocess  the  solution. 

5.5  Save: 


Choosing  the  Save  option  from  the  STORESIM  main  menu  accesses  the  Save  form,  which  is 
used  for  writing  a  copy  of  the  data  base  to  a  user-prescribed  file.  This  file  can  be  used  to  post-process  the 
solution  and/or  to  continue  the  execution  of  the  flow  simulation. 

5.6  Done: 


Chosing  the  Done  option  from  the  STORESIM  main  menu  the  session  can  be  stopped. 

Run  Tetmesh:  Once  the  data  is  stored  in  the  entered/converted  to  ESS  format 

6.0  TETMESH  Command  Line  Options: 

The  valid  options  to  mn  stand  alone  version  of  TETMESH  are.  These  options  are  self 
explanatory.  The  syntax  to  execute  stand  alone  TETMESH  is  :  STORESEVI_DIR/tet  [option]  file.e3s 
[option] 

Valid  options  ( [short-form]  [long-form] ) ; 

[-A]  [ — allow-epsilon-boundary-node-movement] 

[-a  x]  [ — aspect-ratio-threshold=x] 

[-b]  [ — batch] 

[-B]  [ — brute-force-edge-check] 

[-C  x]  [ — curvature-measure=x] 

[-d]  [ — dump-boundary-points] 

[-e  x]  [ — edge-delta=x] 

[-f  x]  [ — flatness-measure=x] 

[-g]  [ — ginternS-debug] 

[-h]  [ — help] 

[-i  x]  [ — ^max-interior-nodes=x] 

[-1]  [ — interactive] 

[-n  x]  [ — max-nodes=x] 

[-0  x]  [ — output-basename=x] 

[-0  x]  [ — solver-option=x] 

[-r  x]  [ — radius-factor=x] 

[-F  x]  [ — output-format=x] 

[-S  x]  [ — tet-shrink-factor=x] 

[-Q  x]  [ — auto-save=x] 

[-q  x]  [ — auto-save-file=x] 

[-R  x]  [ — restart-f ile=x] 

[-S  x]  [ — restart-step=x] 

[-C  x]  [ — restart-codeID=x] 
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[-t  x]  [ — tolerance-point-same=x] 

[-T]  [ — fp3-rej ection-test-of f ] 

[-V]  [ — verbose] 

[-v]  [ — version] 

[-N]  [ — net-to-e2s] 

7.  Input  Parameters: 

Fluid  properties  and  run  control  parameters  can  be  input  to  the  program  through  an  input  file. 
The  default  fluid  properties  and  odier  job  control  parameters  (solver  option,  time  step  scheme,  restructure 
frequency  etc.,)  are  input  through  the  file  ‘.StoresimDefaults’.  These  values  can  also  be  modified  in  the 
widgets.  All  derived  quantities  such  as  the  Mach  No.,  Reynolds,  No.  cannot  be  modified,  instead  they  are 
recomputed  internally  and  displayed  in  the  widget.  To  make  the  program  read  the  ".StoresimDefaults"  file, 
before  executing  the  job  you  need  to  type:  setenv  COMCODEFAULTSFILE  .  StoresimDefaults 
A  sample  .StoresimDefaults  is  given  below. 


JobName 

Title 

Fluid_Name 

Ref_Density 

Ref_Viscosity 

Ref_Temp 

Ref_Velocity 

Ref_Length 

Gamma 

Gascons tant 

SpecificHeat 

Conductivity 

Sutherland_Const_Tl 

Su  ther 1 and_Cons  t_T2 

cfl 


spheres . e3  s 

Flow  over  two  spheres 

AIR 

1.225 

1.7876E-05 

288.15 

1020.8823 

1.0 

1.4 

287.05307 

1004.675 

0.02624 

110.0 

120.0 

0.3 


nsteps_initial  100 

nsteps_simulation  10 

Solver_Option  2 

Time_Integ_Scheme  1 

Use_Local_Or_Global_Time_Step  0 
Restructure_Freguency  5 

Total_Time_Steps  50 

Adaptation_Frequency  0 

Save_Step_Frequency  25 

Max_Mem_In_A_Step  0 

Delta_T  0.025 

Invoke_Solver_Frec3uency  1 

Pseudo_CFL  0.2 

Ps  eudo_Mach  0 . 4 

Pseudo_Reynolds  0.5 
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Display  the  edges  of  selected  surfaces 

Display  the  names  of  selected  surfaces. 

Display  the  edge  names  of  selected  surfaces. 

Display  the  wireframe  mesh  for  selected  durfaces. 

Plot  the  wireframes  of  selected  surfaces  in  non-uniform 
parameter  space. 

Display  the  edge  discretization. 

Display  surface  triangles/nodes. 

Not  Used. 

Display  unenforced  faces. 

Display  surface  normals. 

Reverse  the  orientations  of  the  surface  normals  of 
selected  surfaces. 

Display  the  boundary  conditions. 

Not  Used 

Display  all  elements  with  aspect  ratio  <  the  value  given 
for  Aspect  (If  all  elements  have  aspect  ratio  >  Aspect, 
ail  elements  get  displayed). 

Display  all  elements. 

Minimum  allowable  aspect  ratio. 
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Select  all  surfaces  in  a  given  component 


Select  all  surfaces  in  the  geometry. 
De-select  all  selections. 


Open  the  surface  spacing  parameter  editor. 

Reset  the  graphics  viewport  to  default  position. 
Open  the  Output/Save  control  form. 


Select  all  edges  in  a  selected  component. 
Select  all  edges  in  the  geometry. 
De-select  all  selections. 

Open  the  edge  editor  widget. 

Exit  from  the  TETMESH. 


Do  Everything  (Surf) 
Discretize  Edges 


Triangulate  Surfaces 


mi 


One  step  option  to  triangulate  the  surfaces. 

Discretize  the  edges.  (To  define  the  node  distribution  on  the 
edges  use  the  Edge  Editor  before  discretizing  the  edges). 

Triangulate  all  surfaces. 

Triangulate  only  selected  surfaces. 

Triangulate  only  boundaries  (edges). 

Clear  the  triangulation. 


Do  diagonal  swapping  to  produce  Delaunay  mesh. 
Does  some  data  structure  setup  for  tetrahedralization. 
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One  step  option  totetrahedralize  all  regions. 

Tetrahedralize  the  surfaces. 

Check  if  the  faces  of  surfaces  match  properly  with  neighbors. 
Not  Used. 

Inset  interior  nodes  and  create  tetrahedral  elements. 

Open  the  Statistics  widget. 


Hove  Conponenbs**«  | 

'  1 

Snoobh  Hesh**«  | 

Flou  Parans*** 

Inib  Floufield 

Fire  Sboreflou  | 

fidoptivitj)  Forn*t,  | 

Open  the  menu  to  run  the  six-dof  solver. 

Open  the  menu  to  set  parameters  for  mesh  smoother. 
Display  the  flow  parameters. 

Initialize  the  flow  to  inflow  conditions. 

Execute  the  flow  solver. 

Not  Used. 
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Save  Using  Fornab 
E3s  Fornab 


Generic 


ior'nat 


Oubpub/Save  Conbrol 


Save  Inbo  Exisbing  Or  Neufile:  your_file.e3s 
Save  Inbo  Heiifilej  your^file.gan 
Save  Inbo  Neufile;  your.file*gan-neubral 
Save  Inbo  Hewfilet  your„file*gen 
Save  Inbo  Neufile;  your.file*up 
Save  Inbo  Neufile:  your_f ile*buf 


COHCO  Resbarb  Fornab 


Done 


Surface  Spacing  Paraneters 


Boundary  Spacing  Factor: 
Interior  Spacing  Factor: 
Ji  Spacing 


j  Use  User  Spacing 


I  d  Curvature 


□  Flatness 


I  User  Specified  Spacing: 


1  Flatness  Measure: 

!•  i 

1 

1  Curvature  Measure: 

0 

1 

1  MIH  refinenent  criterion:  | 

i* 

I  MAX  refinenent  criterion:  i 

0 

1  AVE  refinenent  criterion: 

0 

1  Surface  Spacing  Mode: 

1  Surface  Insertion  Mode: 

Ihseri 

I  debugging_level: 


I 


Reset  To  Defaults 


Apply 


Cancel 


Parameters  to  control  the 
triangle  density  at  the  edges 
and  on  the  surfaces.  (Smaller 
value  produces  more  triangles) 

When  selected  uses  above 
parameters. 


The  triangulation  can  also  be 
done  using  any  combinations  of 
these  methods.  When  using 
these,  input  a  non-zero  number 
(appropriately  chosen)  for  the 
corresponding  measure. 


These  are  READ  ONLY 

Valid  Options  are  UV,  XYZ, 
and  XYZ  2 


Do  not  accept  any  changes 

Accept  the  entered  changes 
Close  the  widget. 
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Apply  *= I 


Spacing  Bebueen  Nodes: 
Rppla 


Edges  To  Edit 


Distribution  ;  Off 


Input 


Unif  orn  Cluster 


Cluster  Type  t  ^  Sine 


Hyp,  Tan 


Double  Sine 


Endpoint  Spacing 


Starting  Endpoint  Spacing  : 


Ending  Endpoint  Spacing  : 


J  Display 


Apply 


EKbruded  Thickness  (Along  Nornals) 


Starting  Hornal  Spacing  :  I 


Ending  Hornal  Spacing 


J  Display 


Apply 


Done 
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sphere 

isBdof 


Conponenbs  To  Enploy 


a 

1  □  Extrude 

■II  ■lllliri  1 

Hunber  Layers: 


I  Hun  Passes 
Relax  Factor 
Init  Hun  Passes 
ilnib  Relax  Factor 


jo 

1 

1 

Conponent  Hanes 


Display  Level 


lo 


Reset  To  Defaults 


Extrude  Conponents 


Create  Prisns 


if 


Component  selector. 

All:  To  select  all  components 
De-sel:  To  de-select  all  components. 
To  select  individual  component,  click 
on  its  name. 

I  Extrude  ON/OFF  button 


5 1  Number  of  layers  of  elements  to  be  put 
in  extruded  region. 

Weighting  factor  for  biasing  the  layers. 


Number  of  passes  and  relax,  factor  for 
I  normal  calculations. 


Relax,  factor  and  number  of  passes  for 
averaging  the  normals. 


I  Not  Used 


'  I  To  display  a  given  layer. 
Do  not  accept  changes. 


’'  "I  I  Action  button  to  extrude  component(s). 


I  1  Action  button  to  create  prisms. 


Apply  \  ^Accept  the  entered  changes. 


To  extrude,  follow  this  sequence:  Select  a  component  —  Exturde  ON  -  Apply  —  Extrude  Components 
To  display  a  layer:  Enter  layer  number  in  Display  level  and  click  on  Apply. 
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Move  Conponenbs 


Smooth  Mesh  Parameters . 


Nunber  Of  Tine  Steps; 
Restucbure  Frequency: 
fldapbabion  Frequency: 
Save  Step  Frequency: 


Maxinun  Henory  Per  Step:  |  0 


Tine  Step  Size: 
Solver  Frequency: 

Apply 


0*025000 


Cancel 


Edit  Hoc  Programs... 

Move  Conponenbs 


Display  the  mesh  smoothing 
parameters  w  idg et. 

Max.  number  of  time  steps 

Frequency  of  time  steps  for  mesh 
restructuring. 

Not  Used 

Frequency  of  time  steps  for 
saving  the  solution. 

Not  Used 

Not  Used 

Time  stpe  frequence  for  invoking 
flow  solver. 

Accept  the  entered  changes. 
Close  the  widget 

Open  the  HOC  editor. 

Action  button  to  execute  the 
solver. 


Elenenb  To  Vieu  1 1 


Elem.  number 


Plob  Elenenb 


Plob  Neighbour  Prinb  Grid  Info 


Plot/Print  Options 


Apply 


Cancel 


Accept  the  input. 
Cloe  the  widget 
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Hesh  Snoothing  Paranebers 


I 


Nunber  of  Snoothing  Passes; 
Snoothing  Blgoribhn; 
Relaxation  Factor; 

Convergence  Criterion; 
Hovenenb  Threshold; 
Snoothing  Tine  Step; 

Pseudo  CFL  Nunber; 

Pseudo  Hach  Nunber; 

Pseudo  Reynolds  Nunber; 

Bppia  \ 


I  0*400000 


i 


0*500000 


0*100000 


Cancel 


n6de_floh a  11 


j  0*000000 

1 

0*000000 

1 

. .  . . ,  .i 

1 1*000000 

1 

1  . . 

. ■•'M 

1  0*200000 

. I 

Move  the  nodes  using  an 
incomp,  flow  analogy 
* 


Used  for  solving  the  node 
movement  using  incomp¬ 
ressible  flow  analogy 


I  Accept  the  entered  changes 
Close  this  widget 


Snooth  Mesh 


. i , 


*  Used  for  stand  alone  version  of  mesh  smoother  only. 
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Network  File  Format 


A 
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§ 
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0.984808  0.173648  0 
0.939693  0.34202  0 
0.866025  0.5  0 
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Ed^e  Extrusion  (E3S  Format 
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Boundary  Condition  Specification 
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/*  fai— field  */ 
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Test  Case  No.  1: 


Name:  Concentric  Spheres 
Description; 

This  case  consists  of  a  spherical  rigid  body  of  unit  diameter  embedded  in  a  spherical 
flow  domain.  The  diameter  of  the  outer  sphere  is  10  units.  The  purpose  of  this  test  case  is 
twofold; 

1.  Provide  a  basic  exercise  in  mesh  generation  and  flow  boundary  condition 
specification  using  TETMESH. 

2.  Create  an  integrated  rigid— body  dynamics  and  fluid  flow  simulation  example. 

Features:  General,  Multiple  bodies.  Coupled  fluid  flow,  rigid  body  motion,  and  mesh 
restructuring. 

Input  File:  ./sphere/sphere-patch.net 
Mesh  Generation  Procedure  Outline: 

a.  Convert  the  network  file  to  E2S  format.  (Resulting  file:  net2e2s) 

b.  Convert  net2e2s  file  to  E3S  format. 

c.  Impose  boundary  conditions.  (Edit  the  E3S  file).- 

d.  Run  TETMESH  to  generate  the  mesh.  (See  Flowchart). 
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Test  Case  No.  2: 


Name:  Flow  over  a  Delta  wing 

Description:  This  case  consists  of  a  Delta  wing  attached  to  a  wall  at  the  root.  The  purpose 
of  this  test  case  is  to  mesh  the  flow  domain  using  TETMESH  and  analyze  the  flow  field 
using  STOREFLOW.  The  incident  flow  in  this  case  is  at  Mach  2  and  is  at  angle  of  10 
degrees. 

Features:  General,  trimmed  patches,  flow  solver  parameter  specification,  typical  flow 
simulation. 

Input  File:  ./Delta_wing/delta.e3s 
Simulation  Procedure  Outline: 

a.  Familiarize  with  the  boundary  conditions  format.  (Read  the  E3S  file). 

b.  Run  TETMESH  to  generate  the  mesh.  (See  Flowchart). 

c.  Initialize  the  flow  field  to  Inflow  conditions. 

d.  Run  the  flow  solver.  (Refer  to  TETMESH  widget,  Fire  Storeflow  button) 

e.  Postprocess  the  existing  flow  solution.  (Readin  the  save  files) 


Full  Domain 


i  STORE5IM 

WeaSeo  -1  CC 


s; 


Mesh  Near  me  Wing 


STORE5p4 
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Test  Case  No.  3: 


Name:  Cube  with  an  extended  zero  thickness  surface  inside  a  cube 

Description:  In  this  test  case  a  smaller  cube  is  placed  inside  a  larger  cubical  flow  domain. 
The  inner  cube  has  an  extended  flap  which  is  a  zero  thickness  surface.  The  main  purpose 
of  this  exercise  is  to  demonstrate  the  capability  of  TETMESH  to  recognize  and  appropriatly 
treat  such  surfaces. 

Features:  General,  handling  zero  thickness  surfaces. 

Input  File:  ./zero_thickness/cube-0-in-box.e3s 

Mesh  Generation  Procedure  Outline: 

a.  Familiarize  with  the  specification  of  zero  thickness  surfaces. 

b.  Plot  normals  and  note  the  double  sidedness  of  the  zero  thickness  surface. 

c.  Run  TETMESH  to  generate  the  mesh.  (See  Flowchart). 
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Test  Case  No.  4; 


Name:  Firmed  store 

Description:  This  test  case  has  a  relatively  complex  geometry.  It  comprises  of  several 
multiply  connected  surfaces.  It  provides  a  good  example  for  familiarizing  oneself  with 
several  basic  features  of  STORESIM/TETMESH. 

Features:  General,  Multiple  bodies.  Generation  of  trimmed  patches,  coupled  flow  solver, 
with  risid  body  motion  and  mesh  restructuring. 

Input  File:  ./finned_store/wpfsurf.net 
Mesh  Generation  Procedure  Outline: 

a.  Convert  the  network  file  to  E2S  format.  (Resulting  file:  net2e2s) 

b.  Convert  net2e2s  file  to  E3S  format. 

c.  Impose  boundary  conditions.  (Edit  the  E3S  file). 

d.  Run  TETMESH  to  generate  the  mesh.  (See  Rowchart). 

e.  Postprocess  the  existing  flow  solution.  (Readin  the  save  files) 
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